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Abstract 

This  report  describes  work  that  was  done  under  AFOSR  Contract  Number  FA9550-1 1-1-0056,  studying 
the  structure  of  a  model  urban  boundary  layer  flow.  The  model  geometry  consisted  of  a  set  of  plexiglass 
blocks,  and  the  flow  around  this  geometry  was  studied  both  experimentally  as  well  as  computationally. 
For  the  experiment,  a  Stereoscopic  Particle  Image  Velocimetry  (SPIV)  method  was  developed  that  allows 
for  a  three-dimensional  description  of  this  urban  flow,  and  helps  gain  insight  into  the  characteristic  flow 
structures  in  the  streets  and  canyons  of  our  model  urban  geometry.  On  the  computational  side,  a  new  spectral- 
element  code  was  developed  that  was  demonstrated  to  produce  accurate  results,  and  can  scale  to  thousands 
of  processors  on  large  high-performance  computing  systems.  Good  agreement  between  the  experiment  and 
computation  was  demonstrated. 

Most  notably,  wind  tunnel  experiments  were  performed  at  a  number  of  different  angles  of  incidence, 
providing  for  the  first  time  a  detailed  overview  of  the  effect  of  wind  direction  on  the  flow  structures  in  the 
urban  geometry.  Valuable  information  about  the  flow  structures  are  presented.  The  effects  of  incidence 
angles  from  0  to  45  degrees  of  the  incoming  flow  with  respect  to  the  urban  array  are  investigated.  A  major 
observation  from  this  work  is  that  a  strong  channeling  effect  is  observed  for  all  incidence  angles  and  is  in 
agreement  with  that  observed  in  other  investigations  for  as  little  as  4  degrees.  This  channeling  significantly 
affects  the  turbulence  distribution  within  the  array,  the  correlations  between  the  various  gust  components 
and  the  structures  responsible  for  contaminant  transport. 


Chapter  1 


Introduction 


1.1  Introductory  comments 

Urban  flows  have  been  increasingly  studied  during  the  past  decade.  Many  challenges  are  yet  to  be  fully 
addressed,  with  the  understanding  of  contaminant  dispersion  being  one  of  the  critical  issues.  In  the  context 
of  pollution  dispersion  in  densely  populated  areas  and  the  fear  of  chemical  or  biological  attacks,  more  studies 
are  addressing  this  type  of  flow.  The  physics  involved  are  very  complex  since  the  flow  is  often  strongly 
three-dimensional  and  dependent  upon  the  temperature  distribution,  the  turbulence  induced  by  moving  cars, 
the  presence  of  obstacles  such  as  trees  in  the  streets,  etc.  Remote  flow  sensing,  in  particular,  has  become 
a  very  important  field  of  research.  Through  the  use  of  sparse  sensors  optimally  distributed  in  an  urban 
environment,  the  desire  is  to  estimate  the  characteristics  of  the  flow  field  and  predict  its  evolution.  This  is 
especially  important  in  the  context  of  predicting  the  displacement  of  a  plume  of  contaminants. 

Another  relatively  recent  area  of  interest  relevant  to  this  type  of  study  is  the  desire  to  fly  Micro  Aerial 
Vehicles  (MAVs)  in  urban  areas.  These  MAVs,  which  are  unmanned  aerial  vehicles  with  relatively  small 
dimensions  (10-20  cm  wingspan),  can  be  used  for  surveillance,  data  collection  (remote  sensing),  and  can 
be  flown  in  areas  that  are  hazardous  to  people.  However,  they  are  often  still  limited  to  flying  in  relatively 
calm  wind  conditions  (Watkins  et  al.  [2009,  2010]).  Therefore,  there  is  a  growing  need  to  characterize  wind 
gusts  and  both  the  intensity  and  spatial  distribution  of  turbulence  at  the  urban  scale  so  as  to  provide  the  flow 
fundamentals  to  design  new  generations  of  MAVs  with  enhanced  maneuverability. 

Turbulent  fluid  flows  in  complex  domains  occur  in  nature  and  in  many  industrial  applications.  A  good 
understanding  of  flow  physics  as  well  as  the  ability  to  accurately  predict  these  flows  play  an  important 
role  in  many  applications  ranging  from  weather  prediction  to  efficient  design  of  engineering  equipment. 
Traditionally,  given  the  exorbitant  number  of  grid  points  required  for  accurate  resolution  of  all  flow  features, 
experimental  measurements  supplemented  with  theory  were  the  only  feasible  choice  for  understanding  these 
flows.  Experiments  are  expensive  and  often  provide  data  only  in  a  limited  region  of  the  flow  field.  However 
with  rapid  increase  in  computing  power  and  advances  in  numerical  algorithms,  numerical  simulations  are 
increasingly  becoming  feasible  for  higher  Reynolds  number  flows.  While  resolving  all  scales  of  interest  is 
still  intractable  for  flows  of  engineering  interest,  large  eddy  simulations,  which  resolve  the  large  scales  and 
model  small  scales,  are  becoming  increasingly  feasible. 

Accurate  numerical  simulations  of  turbulent  flows  in  domains  of  engineering  interest  require  algorithms 
which  have  low  dispersion  and  diffusion  errors,  support  complex  geometries  and  are  highly  scalable.  While 
spectral  methods  do  an  excellent  job,  their  application  is  limited  to  simple  geometries.  Second-order  finite 
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volume  and  finite  element  codes  can  handle  complex  geometries  and  are  currently  the  standard  tools  of 
choice  in  the  CFD  industry.  However,  they  have  high  diffusion  and  dispersion  errors  and  need  more  degrees 
of  freedom  for  the  same  accuracy,  making  them  very  expensive  for  long  term  integration.  Spectral  elements 
(Deville  et  al.  [2002],  Karniadakis  and  Sherwin  [2005])  combine  the  high-order  accuracy  of  spectral  methods 
with  the  geometric  flexibility  of  finite  element  methods.  They  have  very  low  dispersion  errors  which  make 
them  very  efficient  for  long  term  integration.  Weak  67°  coupling  between  adjacent  elements  increases  the 
computation  to  communication  ratio  which  facilitates  highly  scalable  algorithms. 

The  computational  work  for  this  project  focused  on  the  development  and  optimization  of  a  parallel 
spectral  element  solver,  Spec  solve,  for  simulating  turbulent  flows  in  complex  geometries.  The  code  uses 
the  Vn-Vn-2  formulation  for  the  spatial  discretization  of  velocity  and  pressure.  Velocity  is  represented 
within  each  spectral  element  using  a  tensor  product  of  Lagrangian  basis  functions  based  on  GLL  (Gauss- 
Lobatto-Legendre)  nodes,  and  C°  continuity  is  enforced  between  adjacent  spectral  elements.  Pressure  within 
each  spectral  element  is  represented  using  a  tensor  product  of  one-dimensional  Legendre  polynomials  of 
appropriate  order  and  inter-element  continuity  is  not  explicitly  enforced.  The  code  is  built  in  C++  and  uses 
the  Message  Passing  Interface  (MPI)  for  communication.  One  of  the  principal  bottlenecks  in  the  numerical 
solution  of  the  Navier-Stokes  equations  is  the  pressure  solve.  A  multilevel  strategy  is  implemented  to  build  a 
scalable  solver  for  pressure.  A  state-of-the-art  algebraic  multigrid  (AMG)  solver  is  implemented  to  solve  the 
coarse  pressure  problem.  An  FDM  (fast  diagonalization  method)  based  block-Jacobi  preconditioner  is  used 
for  efficient  solution  of  the  fine  pressure  problem.  The  code  uses  filter-based  stabilization  for  simulating  high 
Reynolds  number  flows.  A  novel  turbulent  outflow  boundary  condition  is  implemented  in  order  to  stabilize 
flow  at  the  outflow  boundary  in  the  presence  of  localized  regions  of  reverse  flow.  This  parallel  spectral 
element  solver  is  used  for  simulation  of  incompressible  flows  in  model  urban  street  canyons. 

1.2  Advantages  and  limitations  of  field  experiments  and  numerical 
simulations 

While  they  are  essential  to  verify  results  obtained  both  from  numerical  and  wind  tunnel  experiments,  direct 
field  studies  are  extremely  difficult  and  expensive  to  conduct.  Many  parameters  are  influencing  the  flow 
field  and  results  are  often  difficult  to  interpret.  In  terms  of  practicality,  the  most  common  tools  used  to 
measure  velocity  in  field  experiments  are  sodars  and  anemometers  installed  on  tall  towers  which  usually 
results  in  a  rather  coarse  spatial  resolution.  As  a  result,  turbulence  characterization  is  mostly  carried  out  in 
ID,  usually  along  a  vertical  axis,  see  Nielsen  [2000]  and  Christen  et  al.  [2003]  for  example.  Therefore,  no 
direct  information  on  the  spatial  structure  of  the  flow  in  the  directions  normal  to  the  line  of  measurement  is 
available.  Few  field  experiments  have  considered  spatial  variation  of  turbulence  at  the  street  level  (Rotach 
[1995],  Eliasson  et  al.  [2006]).  Other  types  of  work  in  simplified  field  experiments  can  be  found  in  Louka 
et  al.  [1998],  who  worked  toward  a  better  understanding  of  the  ventilation  of  pollution  from  a  street  formed 
of  two  long  buildings  and  its  dependence  on  the  shape  of  the  roofs  of  the  obstacles.  They  focused  on  the 
coupling  of  the  turbulent  airflow  in  a  street  canyon  with  the  turbulent  airflow  above  the  roofs.  Within  the 
same  field  setting,  Louka  et  al.  [2000]  performed  experiments  to  study  the  air  flow  within  the  street  and 
found  that  the  recirculation  region  is  highly  unsteady  and  that  the  shear  layer  from  the  upstream  roof  is  very 
unstable  due  to  a  Kelvin-Helmholtz  instability.  In  a  more  realistic  environment,  Louka  et  al.  [2002]  looked 
at  the  thermal  effects  on  the  flow  field  in  a  street  in  Nantes,  France,  and  compared  their  results  to  a  2D 
numerical  simulation  using  a  standard  k-e  model.  They  found  that  the  numerical  simulations  overestimate 
the  thermal  effects  on  the  flow  field.  Dobre  et  al.  [2005]  showed  the  flow  behavior  in  a  few  streets  in  London, 
UK.  More  specifically,  their  work  illustrated  the  effect  of  the  wind  incidence  angle  with  respect  to  the  streets 
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and  the  role  of  the  street  intersections  on  the  wind  direction  switching  phenomenon.  Their  measurements 
showed  that  the  flow  within  the  street  can  be  viewed  as  the  vector  sum  of  a  channeling  in  the  streets  and  a 
recirculation  vortex. 

The  work  by  Xie  et  al.  [2003]  is  very  interesting  in  the  sense  that  it  would  be  extremely  difficult  to 
reproduce  in  a  wind  tunnel.  They  performed  field  tests  with  the  measurements  of  the  concentration  of  car 
exhaust  in  Guangzhou,  China.  They  found  that  in  the  street,  pollutants  are  carried  from  the  windward  to 
the  leeward  side  of  the  street  by  a  large  recirculation  region.  One  part  of  the  pollutants  gets  trapped  in 
the  recirculation  region,  while  the  other  is  flushed  out  of  the  street  at  the  roof  level.  One  major  finding  is 
that  photochemical  reactions  involving  O3  take  place  at  the  roof  level  and  is  therefore  less  of  a  threat  to 
pedestrians. 

Recently,  there  has  been  an  increased  effort  in  developing  and  improving  computational  and  numerical 
models.  Camelli  et  al.  [2006]  have  used  very  large-eddy  simulation  (VLES)  to  study  the  release  of  a  passive 
scalar  contaminant  in  realistic  urban  areas.  This  technique  is  equivalent  to  standard  large-eddy  simulation 
(LES),  but  the  largest  scales  are  on  the  order  of  a  few  kilometers.  Krajnovic  and  Davidson  [2000]  have  used 
LES  on  a  single  cube  while  Baik  and  Kim  [1999]  have  used  a  2D  k-e  turbulence  model  for  a  2D  street. 

Kim  and  Baik  [2004]  used  a  3D  numerical  simulation,  renormalized  group  k-e  scheme,  to  solve  the  flow 
within  an  array  of  cubes.  They  investigated  the  dependence  of  the  vortical  structures  on  the  wind  incidence 
angle  and  then  classifed  the  flow  in  three  regimes  depending  on  this  incidence  angle. 

More  recently,  Coceal  et  al.  [2006]  performed  a  direct  numerical  simulation  of  the  turbulent  flow  within 
an  array  of  cubes  and  found  very  good  agreement  with  experimental  data.  Only  numerical  studies  can 
provide  information  with  excellent  spatio-temporal  resolution,  but  unless  DNS  is  used,  the  main  drawback 
is  their  limited  ability  to  capture  accurately  the  range  of  turbulence  scales.  The  use  of  DNS  is  however 
computationally  prohibitive.  Even  with  recent  progress  in  computing,  numerical  simulations  need  validation 
through  experimental  studies  for  such  complex  flow  fields.  Therefore,  there  is  a  very  specific  need  for 
reliable  3D  data  in  a  more  realistic  urban  model. 

The  study  of  Kim  and  Baik  [2003]  is  interesting  in  the  sense  that  by  using  a  2D  numerical  model  to  solve 
the  flow  within  a  single  street,  it  allowed  them  to  vary  the  inflow  turbulence  intensity  easily.  Such  a  study 
in  an  experimental  setting  would  be  more  challenging.  They  found  that  as  the  inflow  turbulence  increases, 
both  turbulent  kinetic  energy  and  turbulent  diffusivity  increase.  Similarly,  the  magnitude  of  the  mean  speed 
increases  and  the  vortical  regions  are  strengthened.  An  increase  in  turbulent  intensity  is  directly  linked  to  an 
enhanced  pollutant  dispersion  within  the  street. 

As  mentioned  previously,  thermal  stratification  of  the  atmospheric  boundary  layer  is  very  difficult  to 
simulate  in  a  wind  tunnel  although  some  have  been  successful  in  doing  so  as  will  be  seen  in  the  following 
section.  However,  there  is  no  such  constraint  for  numerical  simulations,  and  researchers  have  started  looking 
into  the  effects  of  solar  radiation  on  the  flow  field  at  the  street  level.  By  simulating  ground-level  heating 
through  various  numerical  models,  Xie  et  al.  [2006]  were  able  to  show  that  thermal  effects  could  enhance 
vertical  diffusion  especially  for  low  mean  wind  speeds.  In  another  work  by  Xie  et  al.  [2005],  the  study  of 
step-up  and  step-down  configurations  for  the  street  geometry  was  carried  out  as  well  as  the  heating  of  either 
the  windward  or  leeward  side  of  the  street.  The  thermal  effects  in  these  cases  could  then  either  strengthen 
or  weaken  the  main  vortex  within  the  street.  Similar  results  were  obtained  independently  by  Mestayer  et  al. 
[1995].  Uehara  et  al.  [2000]  looked  at  the  effect  of  flow  stratification  on  the  flow  field  within  the  street  over 
a  range  of  Richardson  numbers  (defining  stable,  neutral,  or  unstable  atmospheres).  Kim  and  Baik  [2001] 
performed  a  numerical  simulation  of  a  street  with  bottom  heating.  They  were  able  to  observe  significant 
changes  in  the  flow  structure  for  large  temperature  gradients  with  the  generation  of  a  second  vortical  structure 
within  the  street,  located  directly  underneath  the  vortical  structure  found  for  no  temperature  gradients. 

Hang  et  al.  [2012]  also  used  CFD  methods  to  clarify  the  role  of  building  layouts  and  height  variability 
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in  ground-level  pollutant  dispersion.  They  could  distinguish  between  the  roles  of  the  mean  flow  and  the 
turbulent  diffusion  in  pollutant  dispersion  for  each  of  their  building  geometries.  They  also  partially  used 
wind  tunnel  hot-wire  data  to  verify  their  CFD  results. 


1.3  Wind  tunnel  experiments 

Unlike  field  experiments,  wind  tunnel  experiments  allow  for  good  control  over  the  parameters  driving  the 
flow  field  of  interest.  It  also  allows  for  independently  investigating  the  effect  of  the  key  parameters. 

Most  past  wind  tunnel  studies  have  focused  their  interest  on  two-dimensional  (2D)  configurations  (i.e.  span¬ 
ning  the  width  of  the  wind  tunnel)  and  have  considered  the  neutrally  stratified  atmospheric  boundary  layer 
since  it  is  extremely  difficult  to  simulate  a  temperature  distribution  in  a  wind  tunnel.  A  seminal  study  of 
flow  patterns  in  an  urban  environment  was  performed  by  Oke  [1988].  This  work,  which  is  relevant  to  street 
design  and  contaminant  dispersion,  described  the  effect  of  stream  wise  spacing  on  street  canyons.  Although 
his  results  were  2D,  he  identified  three  different  regimes  directly  dependent  on  the  stream  wise  spacing  of 
the  street  canyon.  Li  et  al.  [2008]  performed  experiments  in  a  water  channel,  investigating  the  different  flow 
regimes  using  Laser  Doppler  Anemometry.  Similar  regimes  have  been  found  in  3D  configurations  (i.e.  finite 
size  of  the  streets).  Martinuzzi  and  Havel  [2000]  investigated  the  effect  of  the  stream  wise  spacing  between 
two  cubic  blocks  mounted  in  tandem  in  a  thin  boundary  layer  and  found  three  distinct  regimes  for  differ¬ 
ent  stream  wise  spacings.  Other  parameters  such  as  the  span-to- width  ratio  of  the  obstacles  or  the  span  wise 
spacing  for  obstacle  arrays  are  believed  to  have  an  effect  on  the  transition  between  regimes.  Few  wind  tun¬ 
nel  experiments  have  investigated  more  realistic  geometries  (Rafailidis  [1997],  Kastner-Klein  and  Rotach 
[2004]).  While  these  studies  are  considering  more  complex  geometries,  the  amount  of  information  about  the 
mean  flow  field  and  the  turbulence  within  the  model  is  limited  spatially. 

Recently,  researchers  have  begun  investigating  the  effect  of  three-dimensionality;  but  again,  there  are 
very  few  measurements  offering  full  3D  sets  of  data.  Becker  et  al.  [2002]  with  Laser  Doppler  Anemometry 
(LDA)  and  Sousa  [2002]  with  2D  -  three  component  particle  image  velocimetry  (2D-3C  PIV)  measurements 
were  among  the  first  to  have  looked  at  three-dimensional  flow  fields  around  a  single  cuboid  obstacle.  Becker 
et  al.  [2002]  described  the  effect  of  the  angle  of  attack,  the  aspect  ratio,  the  Reynolds  number  and  the 
boundary-layer  type  on  the  flow  structures  around  a  single  obstacle.  They  found  no  fundamental  changes 
in  the  vortex  structure  for  the  various  boundary  layers  investigated,  and  determined  that  an  increase  in  the 
power-law  coefficient  n  in 


caused  a  reduction  in  the  size  of  the  recirculation  region  downstream  of  the  obstacle.  Here,  U  is  the  time- 
averaged  velocity  at  the  wall-normal  coordinate  z,  and  h  is  a  wall-normal  location  in  the  upper  part  of  the 
inertial  sublayer  at  which  the  velocity  U(h)  is  known.  They  also  found  that  increasing  the  incidence  an¬ 
gle  created  and  subsequently  amplified  a  dislocation  of  one  leg  of  the  arch  vortex  behind  the  single  block 
until  it  attached  to  its  top  for  incidence  angles  larger  than  60°.  Sousa  [2002],  using  conservation  of  mass 
to  extract  the  third  velocity  component  from  2D  PIV  data  around  a  cube,  focused  his  interest  on  the  iden¬ 
tification  and  localization  of  large-scale  vortical  structures.  He  found  that  swirling  strength  and  normalized 
angular  momentum  techniques  were  more  advantageous  for  identifying  coherent  structures  as  compared 
with  a  vorticity -based  method.  Similarly,  Martinuzzi  and  Tropea  [1993],  using  various  flow  visualization 
techniques  as  well  as  static  pressure  measurements  on  the  surface  of  the  obstacle,  showed  the  effect  of  the 
width-to-height  ratio  on  the  three-dimensionality  of  the  flow  past  a  single  3D  obstacle.  They  found  that, 
for  an  obstacle  with  a  width-to-height  ratio  W/H  >  6,  the  flow  downstream  of  the  block  is  largely  2D 
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with  alternating  saddle  and  nodal  points  on  the  windward  face  of  the  obstacle.  They  concluded  that  cellular 
structures  were  created  upstream  of  the  obstacle  and  that  the  flow  was  following  preferred  paths  over  the 
top. 

Hussein  and  Martinuzzi  [1996]  investigated  the  turbulence  dissipation  rate,  the  production,  convection 
and  transport  terms  and  the  balance  of  the  turbulent  kinetic  energy  transport  equation  for  a  three-dimensional 
flow  around  a  cube  mounted  in  a  channel  using  LDA.  They  identified  various  scales  relevant  to  different 
features  of  the  flow  around  the  cube  (e.  g.  wake,  shear  layer,  horseshoe  vortex).  These  different  scales 
characterizing  the  turbulence  production  and  dissipation  are  a  good  indication  that  the  flow  field  within  an 
urban  area  is  a  multi- scale  problem. 

Belcher  and  Coceal  [2001]  and  Coceal  and  Belcher  [2004]  developed  an  urban  canopy  model  for  mean 
winds  in  urban  areas  that  compares  well  with  data  from  wind  tunnel  experiments.  Interestingly,  they  could 
model  the  minimal  distance  within  an  array  of  buildings  that  is  required  to  reach  an  adjustment  to  the 
inhomogeneous  canopies. 

Many  researchers  have  investigated  the  flow  within  an  array  of  obstacles,  for  example,  regular  and  stag¬ 
gered  arrays  of  cubes  (Castro  et  al.  [2006],  Cheng  and  Castro  [2002],  Reynolds  and  Castro  [2008]).  Their 
primary  goal  was  to  understand  the  features  of  urban-like  boundary  layers,  what  influenced  them  and  how  to 
parameterize  them.  With  the  same  objective  in  mind,  MacDonald  et  al.  [2002]  performed  experiments  with 
a  regular  array  of  cubes  and  measured  characteristic  mean  flow  and  turbulence  statistics  for  urban  areas. 

Robins  and  Castro  [1977]  investigated  the  plume  dispersion  in  the  vicinity  of  a  wall  mounted  cube  by 
releasing  propane  for  different  configurations  (porous  cube,  release  from  a  point  source  at  different  locations, 
release  from  a  stack,  etc...)  The  concentration  of  propane  was  measured  downstream  of  the  cube  and  it  was 
found  that  an  effective  source  height  concept  could  be  applied  to  estimate  the  concentration  field  beyond 
two  cube  heights.  A  fundamental  difference  was  found  in  the  behavior  of  the  flow  dispersion  phenomenon 
between  high  and  low  level  release  with  the  authors  underlining  a  possible  effect  of  surface  geometry  on 
dispersion. 

Castro  and  Robins  [1977]  investigated  the  pressure  distribution  on  the  surface  of  a  wall-mounted  cube 
in  a  turbulent  boundary  layer  and  the  flow  field  in  its  wake.  Two  incidence  angles  of  the  incoming  flow 
were  studied.  They  showed  that  an  increase  of  turbulence  intensity  and  shear  would  reduce  the  size  of  the 
recirculation  region  downstream  of  the  cube.  Reattachment  on  the  top  surface  of  the  obstacle  is  observed  for 
the  urban  type  boundary  layer  (as  opposed  to  the  cube  being  placed  in  a  uniform  flow).  This  reattachment 
occurs  for  different  cube  sizes,  but  may  be  intermittent  for  the  larger  obstacles  (as  the  height  approaches  the 
boundary  layer  thickness).  For  a  45°  incidence  angle,  the  cube  placed  in  a  uniform  flow  sheds  vortices  from 
its  top  edges,  similar  to  a  delta  wing,  that  have  a  strong  effect  on  the  axial  velocity  for  a  distance  of  about  six 
times  the  cube  length,  then  this  effect  diminishes  quickly.  When  turbulence  intensity  is  raised  in  the  case  of 
the  urban  type  boundary  layer,  this  effect  tends  to  disappear.  Garbero  et  al.  [2010]  performed  an  experimental 
study  of  pollutant  dispersion  within  a  street  network,  measuring  the  mass  exchange  between  the  streets  and 
the  flow  above  the  roofs  for  various  geometrical  configurations  as  well  as  various  wind  directions. 

While  the  majority  of  investigations  concentrate  on  street  canyons  in  the  middle  row  of  a  full  array  of 
buildings,  Princevac  et  al.  [2010]  used  PIV  measurements  in  a  water  channel  to  focus  on  the  flow  distribution 
for  a  general  array  of  mock  up  buildings,  they  discovered  a  new  flow  feature  called  ‘lateral  channeling’  which 
explains  the  inward  and  outward  mean  flow  distribution  on  pathways  on  the  boundaries  of  the  urban  array 
model. 

An  interesting  experimental  study  at  a  larger  scale  was  carried  out  by  Richards  et  al.  [2001]  where  they 
measured  the  surface  pressure  on  a  6  m  cube  sitting  on  the  ground  in  an  open  country.  They  found  good 
agreement  with  wind  tunnel  data  for  the  windward  face  of  the  cube,  but  the  pressure  distribution  in  the  field 
experiment  was  more  sensitive  to  incoming  velocity  profiles,  turbulence  and  Reynolds  number  on  the  other 
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faces  than  the  pressure  distribution  obtained  in  the  wind  tunnel  experiment.  Richards  et  al.  [2007]  carried  out 
a  1 :40  scaled-down  wind  tunnel  experiment  of  this  cube  and  compared  their  results  to  the  field  experiments. 
They  found  similar  behavior  but  identified  various  sources  of  discrepancies  between  the  wind  tunnel  and 
field  experiments.  These  included  a  Reynolds  number  effect  and  the  difficulty  of  reproducing  varying  wind 
directions  to  match  the  field  experiment,  affecting  the  mean  and  amplitude  of  pressure  observed  during  the 
measurements. 

Takimoto  et  al.  [201 1]  applied  PIV  in  an  outdoor  field  experiment  on  an  array  of  1.5  m- sized  cuboid  ob¬ 
stacles.  In  spite  of  several  differences  compared  to  the  wind  tunnel  test,  they  suggested  valuable  information 
regarding  turbulent  flow  structures  by  proposing  dominant  eddy  modes  and  the  shear  layer  behavior  at  the 
canopy  roof  level.  A  thorough  study  about  the  effect  of  the  array  geometry  parameters  such  as  array  layout, 
buildings  heights,  and  wind  direction  on  aerodynamic  parameters  of  buildings,  such  as  drag  coefficient  (Cd) 
was  carried  out  by  Hagishima  et  al.  [2009]  in  the  wind  tunnel  using  floating  pressure  gauges. 

Kellnerova  et  al.  [2012]  used  PIV  to  investigate  the  turbulent  flow  in  street  canyons.  POD  methods 
were  used  to  analyze  the  velocity  field  and  extract  coherent  modes.  They  studied  the  effect  of  roof  top 
shape  (pitched  and  flat)  on  turbulent  flow  characteristics  in  streets.  Even  with  valuable  improvements  in 
urban  flow  studies,  extracting  three-dimensional  shape  and  distribution  of  the  vortical  structures  was  still  a 
weak  aspect  of  most  investigations.  Monnier  et  al.  [2010]  for  the  first  time  used  stereoscopic  PIV  (SPIV) 
methods  to  generate  fully  three-dimensional  velocity  data  of  street  canyons  in  wind  tunnel  tests.  They 
used  a  scaled  model  of  the  MUST  field  experiment  (Biltoft  [2001]),  investigated  the  effects  of  streamwise 
spacing,  ambient  boundary  layer  regimes  (wake  interference  and  skimming  flow),  and  small  wind  direction 
angles  (0°,  4.5°  and  15°).  Representation  of  mean  streamlines,  turbulent  kinetic  energy  (TKE)  contours  and 
coherent  structures  were  used  in  the  investigation.  Using  proper  orthogonal  decomposition  (POD)  methods 
to  refine  statistical  data  provided  them  with  advanced  three-dimensional  iso-surfaces  of  coherent  structures 
of  vortical  motions.  Arch  vortices  were  clearly  detected  and  depicted  in  streets.  The  effect  of  wind  direction 
on  the  flow  characteristics,  dispersion,  TKE  and  the  mean  flow  was  studied  by  Kim  and  Baik  [2004]  using 
CFD.  Angles  of  0°  to  45°  in  increment  of  5°  were  studied.  The  evolution  of  the  portal  vortex  (arch  vortex) 
was  studied  in  three  main  categories  of  wind  direction  angle  of  0°,  5°  to  20°  and  25°  to  45°. 

Mostly  during  the  last  decade,  researchers  have  begun  to  investigate  other  effects  on  the  flow  field  within 
urban-like  environments.  The  addition  of  thermal  effects,  presence  of  obstacles  such  as  trees,  and  extra 
turbulence  generated  by  moving  cars  are  more  and  more  investigated  experimentally.  In  a  first  study,  Kovar- 
Panskus  et  al.  [2002b]  looked  at  the  streamwise  aspect  ratio  effect  of  the  street  on  the  flow  pattern,  followed 
by  a  study  where  Kovar-Panskus  et  al.  [2002a]  investigated  thermal  effects  in  a  wind  tunnel  experiment. 
They  investigated  the  flow  field  within  a  street  for  conditions  where  a  single  recirculation  region  is  the 
dominant  flow  structure.  Under  low  speed  wind  conditions,  as  the  heating  of  the  windward  wall  is  increased, 
the  reverse  flow  speed  magnitude  is  reduced.  Upon  further  heating,  the  flow  pattern  significantly  changed. 

Some  researchers  have  started  to  look  into  added  turbulence  generated  by  moving  cars.  For  example, 
Eskridge  and  Rao  [1986]  focused  specifically  on  the  turbulence  levels  created  by  moving  objects.  More 
recently,  Kastner- Klein  et  al.  [2001]  studied  the  high  concentration  levels  of  car  exhaust  encountered  when 
the  incoming  flow  was  normal  to  the  street,  and  looked  at  the  effects  of  one-way  and  two-way  traffic  on 
this  concentration.  In  another  study,  Kastner-Klein  et  al.  [2000]  found  that  dispersion  of  pollutants  could 
be  enhanced  by  the  motion  of  vehicles,  especially  at  low  wind  speeds.  Ahmad  et  al.  [2002]  implemented 
multi-lane  traffic  for  simulating  road  traffic  and  concluded  that  pollutant  concentrations  could  be  reduced 
significantly  by  moving  cars,  once  again  at  low  wind  speeds. 

In  the  last  few  years,  the  effect  of  still  obstacles  on  the  dispersion  of  pollutants  has  been  investigated. 
More  specifically,  the  effect  of  trees  was  highlighted  in  the  work  by  Gromke  and  Ruck  [2009],  Gromke  and 
Ruck  [2007]  and  Gromke  et  al.  [2008].  It  was  found  that  in  some  cases  the  presence  of  trees  lowers  the 
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dispersion  of  pollutants  from  the  leeward  side  of  the  street,  by  weakening  the  large  recirculation  region.  On 
the  windward  side  of  the  street,  pollutant  concentrations  were  less  than  in  the  baseline  case  (with  no  trees), 
due  to  a  larger  amount  of  clean  air  moving  into  the  street  from  the  top.  Gayev  and  Savory  [1999]  focused  on 
trying  to  model  the  roughness  associated  with  stationary  obstacles  in  a  single  street-canyon  model. 
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Chapter  2 


Wind  Tlinnel,  Apparatus  and  Data 
Processing 

2.1  Modelling  of  the  Atmospheric  Boundary  Layer 

2.1.1  Definition 

The  first  step  in  this  investigation  was  to  experimentally  model  the  atmospheric  boundary  layer  (ABL)  for  the 
conditions  of  interest  Stull  [1988].  The  wind  velocity  profile  of  the  ABL  is  often  characterized  by  a  log-law 
in  the  inertial  sublayer.  However,  using  a  power  law,  see  Equation  1.1,  provides  a  simpler  characterization 
of  a  neutrally  stratified  boundary  layer  Plate  [1971]  since  only  one  parameter  (the  exponent)  is  needed  to 
describe  the  wind  velocity  profile  in  the  inertial  sublayer.  Typical  exponents  for  various  types  of  terrains 
can  be  found  in  the  work  by  Davenport  [1965],  and  range  from  0.4  for  large  high-rises  to  0.16  over  a  flat 
surface  such  as  a  lake.  An  ABL  with  a  power-law  exponent  of  0.17,  corresponds  to  flow  over  lower  suburban 
buildings  and  was  used  for  this  investigation. 

2.1.2  Wind-Tunnel  Modeling 

The  boundary  layers  are  modeled  in  a  closed-loop  wind  tunnel  at  IIT,  based  on  the  work  done  by  Nagib  et  al. 
[1974]  and  Gunnarsson  [1974].  The  experiment  is  carried  out  in  the  low  speed  test  section  that  is  capable  of 
mean  free-stream  speeds,  Uo,  up  to  8  m/s.  The  test  section  is  1.2  m  in  span  by  1.7  m  in  height  and  0.635  m  in 
length.  For  larger  test  sections  such  as  this  one,  it  is  very  important  to  consider  the  span  wise  uniformity  of  the 
incoming  boundary  layer.  Nagib  et  al.  [1974]  and  Gunnarsson  [1974]  used  roughness  elements  and  a  counter 
jet  upstream  of  the  test  section  to  tune  the  characteristic  parameters  of  the  approaching  boundary  layer  and 
ensure  spanwise  uniformity.  The  counter  jet,  discussed  both  in  Nagib  et  al.  [1974]  and  Gunnarsson  [1974], 
proved  to  be  very  useful  in  ensuring  spanwise  uniformity  for  the  free-stream  flow  speeds  of  interest  here. 
The  configurations  we  chose  for  roughness  elements  and  counter  jet  settings  provide  a  spanwise  uniform 
neutrally  stratified  atmospheric  boundary  layer  for  free-stream  speeds  ranging  from  Uo  =  2.2  m/s  to  3.4 
m/s. 

The  counter  jet  consists  of  a  60-mm  diameter  steel  tube  placed  on  the  floor  of  the  wind  tunnel  and 
spanning  the  entire  width.  The  counter  jet  is  placed  at  the  upstream  location  of  the  roughness  fetch  (Figure 
2.1)  which  is  3.48  m  upstream  of  the  test  section.  There  are  38,  6.35-mm  diameter  holes  drilled  along  the 
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span  of  the  steel  tube,  with  the  tube  itself  being  connected  to  a  compressed  air  supply.  The  orientation,  Oj , 
of  the  38  jets  can  be  varied  from  —20°  to  +20°  with  respect  to  the  oncoming  free  stream  by  rotating  the 
tube.  The  magnitude  of  the  jet  velocity,  Uj,  can  also  be  varied  by  controlling  the  compressed  air  supply 
pressure.  For  our  study,  we  found  that  Uj/Uo  =  15  with  an  upstream  flow  angle,  Oj  of  +10°  (see  Fig.  2.1 
for  definition  of  Oj)  provides  a  span  wise  uniform  boundary  layer  in  the  test  section.  The  reader  is  referred  to 
Gunnarsson  [1974]  for  a  more  comprehensive  discussion  of  the  counter  jet  technique.  Directly  downstream 


Counter  Jet 

777777777777 


Roughness  Elements 

Figure  2.1:  Photograph  of  counter  jet  and  roughness  elements  upstream  of  the  test  section.  Note:  the  test 
section  is  not  shown  here  and  is  downstream  of  this  image. 


of  the  counter  jet  is  the  roughness  element  fetch.  The  roughness  elements  used  are  30  mm  cubes  placed 
randomly  on  the  floor  (see  Figure  2.1).  The  end  of  the  roughness  fetch  is  2.08  m  from  the  counter  jet  and 
1.40  m  from  the  upstream  edge  of  the  1.22  m  x  1.93  m  test  section.  The  roughness  elements  occupy  8%  of 
the  planform  area  of  the  entire  roughness  fetch.  The  test  section,  which  begins  1.40  m  after  the  downstream 
end  of  the  roughness  fetch,  can  be  seen  in  Figure  2.2. 

2.1.3  Characteristics  of  the  modelled  ABL 

The  approach  boundary  layer  is  documented  using  a  spanwise  array  of  three  hot  wires  and  its  characteristics 
are  summarized  in  Table  2.1..  The  rack  is  mounted  on  a  vertical  traverse  system  enabling  measurement 
of  the  velocity  profiles  starting  from  a  position  in  close  proximity  to  the  floor  («  1  mm)  and  extending 
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Figure  2.2:  Schematic  of  the  stereoscopic  PIV  setup,  the  roughness  fetch  and  the  test  section. 


Table  2.1:  Boundary  layer  test  matrix 


ABL# 

U0 

s 

e 

Re<9  n 

ABL 

3.4  m/s 

«  450  mm 

«  50  mm 

«  10800  0.17 
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approximately  400  mm  above  it.  The  middle  hot  wire  is  located  along  the  centerline  of  the  wind  tunnel  while 
the  two  other  hot  wires  are  located  127  mm  apart  from  the  centerline.  Note  that  for  meteorological  studies,  x, 
y  and  z  are  the  streamwise,  spanwise  and  wall-normal  coordinates,  respectively.  Since  the  hot-wire  traverse 
is  limited  to  400  mm,  a  Pitot-static  tube  traverse  is  used  to  estimate  the  boundary-layer  thickness  S,  where  S 
is  defined  as  the  wall-normal  (z)  location  where  U/Uo  =  0.99.  The  momentum  thickness  0  is  obtained  from 
the  profiles  (0  «  50  mm)  and  is  computed  as: 


0  = 


U(z) 

U0 


(2.1) 


where  the  integration  is  performed  over  the  interval  0  to  0.8  m. 

Figure  2.3  presents  the  velocity  profile  for  our  ABL.  The  velocity  profile  was  obtained  100  mm  upstream 
of  the  test  section  and  the  values  for  S ,  0  and  the  power-law  coefficient  n  are  presented  in  Table  2.1. 


Figure  2.3:  Logarithmic  representation  of  normalized  velocity  profiles  for  ABL  upwind  of  the  beginning  of 
the  urban  array. 
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2.2  Urban  model 


The  canopy  array  used  in  this  study  consists  of  7  rows  of  blocks  in  span  and  5  in  the  streamwise  direction, 
which  makes  6  parallel  streets  in  span  and  4  intersecting  streets  in  the  streamwise  direction.  The  data  were 
taken  in  the  region  behind  the  mid-span  column  of  the  blocks,  so  the  data  were  acquired  in  the  four  streets 
shown  in  Figure  2.4.  The  three  intersecting  streets  on  each  side  in  span  were  aimed  to  provide  geometrical 
symmetry  and  spanwise  flow  uniformity.  The  geometrical  information  of  the  urban  model  is  shown  in 
Figure  2.5.  The  Cartesian  coordinate  system  and  cardinal  directions  are  used  to  better  address  the  locations 
and  directions  in  the  array.  Note  that  the  coordinate  system  is  attached  to  the  urban  array  in  cases  of  various 
wind  directions.  The  model  dimensions  are  described  as:  building  height,  H  =  50  mm,  building  width, 
W  =  L  =  25  mm,  and  street  length,  S  =  37.5  mm. 


Figure  2.4:  Array  structure  (5x7)  and  data  acquisition  regions,  top  view. 


The  other  parameters  which  describe  the  array  geometry  are  packing  density  and  dimensionless  frontal 
area.  Packing  density  (Ap)  is  the  ratio  of  the  total  plan  area  of  the  buildings  (Ap)  to  the  total  underlying 
surface  area  (At).  Dimensionless  frontal  area  (A/)  is  defined  as  the  ratio  of  the  frontal  area  of  the  buildings 
facing  the  wind  (Af)  over  At.  Dimensionless  frontal  area  (A/)  is  important  in  the  drag  calculations,  since 
it  represents  the  surface  facing  the  wind  flow.  Typical  ranges  of  A  f,  are  from  0.1  for  areas  with  a  moderate 
density  of  buildings  up  to  0.3  in  downtown  areas.  For  the  current  investigation  A/  =  0.4  and  Ap  =  19.9 
%  which  is  considered  a  relatively  medium  dense  array  compared  to  typical  urban  array  models.  Zaki 
et  al.  [2011]  used  packing  densities  of  7.7%  up  to  39%  for  their  experiments,  and  Pillai  et  al.  (2012)  used 
packing  densities  of  6%  up  to  25%,  where  they  were  studying  the  effect  of  the  packing  density  on  pollutant 
dispersion. 

The  desired  region  of  study  in  each  street,  shown  in  Figure  2.6(a),  is  a  volume  of38  x62.5x71  mm,  in 
the  x  x  y  x  z  directions  respectively.  This  region  is  centered  at  the  middle  of  each  street  in  span  and  covers 
half-length  of  intersecting  streets  on  each  side.  In  order  to  attain  a  full  data  matrix  of  this  region,  a  set  of 
data  planes  is  accumulated.  Each  data  plane  is  a  volume  of  38  x  2  x  71  mm, in  xxy  x  z  respectively  (Figure 
2.6(a)).  Each  plane  contains  a  38  x  71  grid  of  data  points.  27  data  planes  are  in  span  with  one  located  at  the 
center  (as  shown  in  Figure  2.6(b))  and  13  on  each  side.  The  spanwise  distances  of  these  13  planes  from  the 
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Figure  2.5:  Geometrical  information  of  the  urban  model. 
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center  on  each  side  are  shown  in  Table  2.2. 


Table  2.2:  Span  wise  data  acquisition  positions. 


Plane  # 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

Distance 

from 

2.5 

4.5 

6.5 

8.5 

10.5 

12.5 

14.5 

16.5 

18.5 

20 

23.5 

27.5 

31.25 

center 

(mm) 

(a)  Region  of  interest  in  a  single  street 


(b)  Single  SPIV  plane 


Figure  2.6:  Data  acquisition  regions. 


In  each  acquisition  plane,  a  set  of  1218  pairs  of  images  are  taken  by  each  camera.  The  camera  shooting 
exposure  time  is  set  to  be  182  /is.  Laser  shooting  frequency  and  camera  capturing  frequency  are  both  15  Hz. 
The  detailed  SPIV  specifications  are  summarized  in  Table  2.3. 


2.3  Hot-wire  anemometry 

The  hot-wire  probes  used  to  characterize  the  three  atmospheric  boundary  layers  described  earlier  were  op¬ 
erated  by  a  Constant  Temperature  Anemometer,  CTA,  custom-made  at  IIT.  The  hot  wires  were  built  using 
3.8-/im  diameter  tungsten  wire  (Sigmund  Cohn  Co.).  The  sensing  length  of  the  wire  was  about  1  mm  to 
ensure  a  length-to-diameter  ratio  of  l/d  >  200.  The  non-dimensional  sensing  length,  l+  =  luT/v ,  is  less 
than  20  as  advised  by  Blackwelder  and  Haritonidis  [1983]  to  satisfy  the  spatial  resolution  of  the  small-scale 
turbulence  structures  found  in  wall-bounded  flows.  The  overheat  ratio  was  set  at  1.7.  The  hot  wires  were 
calibrated  in  situ  against  a  Pitot  static  probe.  Fourth-order  polynomial  curves  were  used  to  fit  the  calibration 
data.  The  temperature  was  also  monitored  using  a  thermocouple  so  as  to  correct  for  any  drift  in  the  measured 
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Table  2.3:  Span  wise  data  acquisition  positions. 


Feature 

Specification 

Laser  type 

pulsed  dual-head  Nd-Yag 

Pulse  separation 

130  /is 

Q-switch 

180  /is 

Laser  frequency 

15  Hz 

Camera  resolution 

1280  x  1024  pixels 

Interrogation  area 

32  x  32  pixels 

Illumination  thickness 

2  mm 

Camera  capturing  rate 

15  fps 

Camera  exposure  time 

182  fi s 

Image  contrast 

0.05 

Number  of  image  pairs/  data  set 

1218 

Number  of  data  sets/  street 

27 

signals  of  the  hot  wires  according  to  the  temperature  correction  formula  provided  by  Drubka  et  al.  [1977]. 
The  output  signals  were  low-pass  filtered  using  an  Ithaco  filter  with  a  2  kHz  cut-off  frequency.  The  sampling 
frequency  was  set  at  4  kHz  and  the  acquisition  time  was  set  at  30  s  per  wall-normal  position  of  the  hot  wire. 
A  National  Instruments  data  acquisition  board  (PCI-6251)  was  used  to  acquire  the  hot-wires  measurements 
along  with  the  Lab  VIEW  interface. 

2.4  Stereoscopic  PIV  system 

Particle  Image  Velocimetry  (PIV)  is  a  measurement  technique  that  was  introduced  in  the  early  80’s,  see 
Adrian  and  Yao  [1985].  Adrian  [2005]  defines  PIV  as  “the  accurate,  quantitative  measurement  of  fluid 
velocity  vectors  at  a  very  large  number  of  points  simultaneously.”  An  extension  to  this  technique  called 
Stereoscopic  PIV  (SPIV)  introduced  by  Soloff  et  al.  [1997]  is  used  in  this  study.  Essentially,  the  addition  of  a 
second  camera  allows  for  determining  the  third  component  of  the  velocity  field  in  the  out-of-plane  direction. 
Figure  2.7  presents  a  simple  sketch  of  the  two  cameras  pointing  at  the  region  of  interest  from  different 
positions.  The  displacement  of  a  particle  is  therefore  seen  differently  from  the  two  cameras.  Through 
a  calibration  procedure  that  records  known  out-of-plane  displacements  of  a  calibration  target  from  both 
cameras,  the  3D  velocity  field  can  be  estimated  in  the  laser  plane.  The  Scheimpflug  condition  (image  plane, 
lens  plane  and  object  plane  intersecting  in  one  line,  see  Louhichi  et  al.  [2006])  is  also  satisfied  so  as  to 
improve  the  image  quality  when  looking  at  the  laser  plane  at  an  angle. 

Figure  2.8  presents  a  photograph  of  the  SPIV  setup  with  the  laser  firing.  In  order  to  make  SPIV  measure¬ 
ments  possible  over  a  large  domain,  it  is  necessary  to  move  the  measurement  system  quickly  and  accurately. 
The  calibration  process  for  SPIV  is  tedious  because  a  new  calibration  is  usually  required  for  each  indepen¬ 
dent  plane  of  data  if  any  of  the  SPIV  components  are  moved.  However,  in  this  work,  a  solution  to  overcome 
this  limitation  is  implemented,  where  we  set  the  entire  SPIV  system  on  a  single  plate  sitting  on  a  two-axis 
traverse  system  located  under  the  wind  tunnel.  As  can  be  seen  in  Figures  2.2  and  2.8,  the  laser  head,  the 
light  sheet  optics,  and  the  two  cameras  are  secured  to  that  plate.  The  stereo  configuration  chosen  is  also 
shown  in  Figures  2.2  and  2.8,  with  the  two  cameras  viewing  the  laser  light  sheet  from  opposite  span  wise 
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True  3D  displacement 


2D  Displacement  seen  from  left 


Right  camera 


2D  Displacement  seen  from  right 


Focal  plane  (center  of  light  sheet) 


Figure  2.7:  Stereoscopic  principle. 


Figure  2.8:  Picture  of  the  Stereoscopic  Particle  Image  Velocimetry  system. 


16 


sides  at  an  angle  /3  of  about  50°  as  defined  in  Figure  2.7.  The  light  sheet  enters  the  test  section  through 
a  thin  (3.18  mm)  glass  section  of  the  floor.  The  advantages  of  this  setup  include  the  high  accuracy  of  the 
computer-controlled  two-axis  traverse  system,  the  ease  of  rotation  of  the  entire  PIV  system  to  study  different 
approach  flow  incidence  angles,  and  most  importantly,  this  system  requires  a  single  calibration  per  dataset 
(up  to  224  planes).  In  addition,  the  calibration  process  is  greatly  simplified  by  keeping  the  calibration  target 
immobile  and  instead  displacing  the  SPIV  setup  at  various  out-of-plane  positions  using  the  traverse  system. 

The  traverse  system  is  a  two-axis  traverse  system  manufactured  by  Velmex,  and  has  a  large  traveling 
distance  (635  mm  along  both  axes)  so  that  the  entire  test  section  can  be  covered  if  needed.  The  repeatability 
is  typically  4  /im,  which  means  that  when  the  laser  plane  is  returned  to  its  initial  position,  there  is  no  more 
than  4  /im  of  discrepancy  due  to  backlash.  The  straight  line  accuracy  is  0.076  mm  over  the  entire  travel 
distance  and  the  screw  lead  accuracy  is  0.076  mm  per  250  mm  (data  provided  by  Velmex). 

The  SPIV  system  is  from  Integrated  Design  Tools,  Inc  (IDT);  the  laser  is  a  pulsed  dual-head  Nd-Yag 
New  Wave  Research  (200  mJ  per  pulse).  The  two  cameras  are  X-Stream  5  with  1280x1024  pixels  and  IDT 
ProVision-XS  software  is  used  to  acquire  and  process  the  raw  images.  The  frequency  of  acquisition  is  always 
set  at  its  maximum  value  of  15  Hz. 

The  seeding  is  done  via  three  Trust  Science  Innovation  (TSI)  atomizers  and  consists  of  a  mixture  of 
polyethylene  glycol  (PEG)  and  distilled  water.  The  atomizers  produce  a  mean  droplet  diameter  of  0.3  fi m 
with  a  geometric  standard  deviation  of  less  than  2.0  /im  so  that  the  particles  will  follow  the  flow  accurately 
(see  Raffel  et  al.  [1998]).  The  seeding  particles  are  injected  into  the  wind  tunnel  through  three  inlets  at 
the  floor  just  downstream  of  the  roughness  fetch.  The  injection  speed  is  low  enough  so  that  it  does  not 
generate  additional  disturbances  in  the  flow.  This  was  confirmed  by  acquiring  hot-wire  velocity  profiles  just 
upstream  of  the  urban  array.  The  atomizers  were  emptied  from  the  water  and  PEG  mixture  and  run  at  the 
same  compressed  air  pressure  value  as  that  used  in  the  SPIV  study  to  simulate  the  seeding.  These  velocity 
profiles  agreed  very  well  with  those  for  no  injection;  therefore,  the  effect  of  the  seeding  on  the  boundary 
layer  is  considered  to  be  negligible. 


2.5  SPIV  data  processing 

The  commercial  software  used  to  process  the  raw  SPIV  images  is  ProVision-XS  by  Integrated  Design  Tools 
(IDT)  and  lets  us  store  a  “flag  matrix”  that  contains  the  status  of  each  computed  velocity  vector.  Namely,  we 
know  if  the  computation  of  the  vector  is  successful;  which  means,  a  valid  vector  is  one  where  all  conditions 
to  obtain  an  accurate  estimation  of  the  velocity  are  met.  Conversely,  a  spurious  vector  is  one  where  all 
conditions  are  not  met.  A  common  approach  is  to  replace  the  spurious  vectors  with  interpolated  vectors 
using  a  least  square  estimation  based  on  the  nearest  neighbors  approximation.  This  is  not  the  approach  used 
for  the  current  study.  When  considering  turbulent  flows  such  as  the  urban  flow,  it  is  very  difficult  to  obtain 
a  uniform  seeding.  Often,  PIV  images  will  present  regions  without  any  particles,  and  depending  on  the  size 
of  these  regions,  the  nearest  neighbors  approximation  can  be  highly  inaccurate. 

A  more  efficient  way  of  recovering  the  missing  information  is  to  use  gappy-POD  (Gunes  et  al.  [2006]), 
which  was  successfully  implemented  by  Murray  and  Ukeiley  [2007]  on  PIV  data.  The  amount  of  spurious 
velocity  vectors  in  our  data  sets  is  on  average  2%  and  the  large  amount  of  SPIV  snapshots  available  to  us  for 
each  vertical  plane  make  gappy-POD  an  excellent  alternative  to  standard  interpolation  Gunes  et  al.  [2006]. 
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2.6  Notes  on  accuracy  and  calibration  of  the  SPIV  setup 

2.6.1  Accuracy  of  out-of-plane  velocity  component 

For  typical  SPIV  investigations,  the  accuracy  for  the  in-plane  components  is  usually  given  as  1-2%.  Recently, 
many  have  looked  into  various  techniques  to  improve  the  accuracy  of  the  out-of-plane  velocity  component 
Calluaud  and  David  [2004],  Lecerf  et  al.  [1999],  ?. 

Prior  to  the  collection  of  the  SPIV  data  in  the  urban  environment,  a  separate  experiment  was  carried  out  to 
estimate  the  accuracy  of  our  SPIV  system.  In  order  to  estimate  the  accuracy  of  the  out-of-plane  component  of 
the  velocity  field,  SPIV  measurements  were  taken  in  an  axisymmetric  jet.  Figure  2.9  presents  a  photograph 
of  the  setup.  The  jet  is  shown  on  the  right  hand  side  in  the  picture  and  the  cameras  are  mounted  on  the  side 
and  pointed  at  a  calibration  target,  which  is  mounted  perpendicular  to  the  jet  axis.  In  this  configuration,  the 
out-of-plane  component  captures  the  stream  wise  velocity  of  the  jet.  Several  data  sets  were  collected  in  this 
configuration  by  varying  the  calibration  target  type  and  the  number  of  calibration  images.  The  calibration 
target  providing  the  best  results  for  this  jet  flow  is  shown  in  Figure  2.10,  as  seen  from  both  cameras.  The 
target  is  made  of  aluminum  to  ensure  its  flatness.  The  grid  nodes  (with  a  5 -mm  spacing  in  both  the  horizontal 
and  vertical  directions)  are  holes  through  the  plate.  In  this  example,  the  white  mesh  displayed  on  top  of  both 
views  of  the  target  is  covering  a  50  x  35-mm  region.  Each  mesh  node  is  aligned  with  a  target  grid  node.  The 
calibration  process  tracks  the  displacement  of  these  grid  nodes  as  the  target  is  moved  by  a  small  amount  in 
the  out-of-plane  direction.  The  best  results  in  accuracy  of  the  mean  velocity  measurements  were  obtained  for 
1 1  calibration  images  spanning  the  thickness  of  the  laser  light  sheet  («  2  mm).  These  calibration  parameters 
were  therefore  used  in  the  collection  of  the  SPIV  measurements  in  the  urban  environment. 

The  same  velocity  measurements  were  performed  for  a  configuration  where  the  calibration  target  was 
mounted  parallel  to  the  jet.  As  a  result,  the  stream  wise  component  of  the  jet  was  captured  by  an  in-plane 
component  of  the  SPIV  data.  Figure  2.11  illustrates  the  results  for  the  stream  wise  velocity  profile  of  the 
axisymmetric  jet  as  measured  using  both  configurations.  It  can  be  seen  that  the  agreement  in  the  mean 
velocity  profile  is  very  good  in  the  core  of  the  jet  (—0.4  <  r/D  <  0.4).  In  this  region,  the  difference  in 
magnitude  between  the  in-  and  out-of-plane  configurations  does  not  exceed  2%  as  shown  in  Figure  2.12  and 
is  comparable  to  the  accuracy  obtained  by  Lecerf  et  al.  [1999].  Larger  differences  are  observed  in  the  shear 
layer  where  the  seeding  was  non-uniform. 

Since  several  manual  steps  are  needed  for  the  calibration  process  in  the  urban  array  study,  seven  sets 
of  calibrations  are  collected  and  compared  for  each  setup,  in  order  to  reduce  the  human  errors.  The  most 
accurate  calibration  should  meet  two  criteria  on  the  out-of-plane  (y  direction)  component  of  the  velocity 
field.  The  candidate  should  show  a  relatively  low  level  of  variance,  and  also  reveal  a  good  agreement 
with  other  calibrations  at  the  same  time.  These  different  calibrations  were  applied  to  street  data  in  order 
to  compare  variance  levels.  Root  mean  squares  (rms)  of  span  wise  velocities  (vrrns)  were  compared  on  a 
vertical  line  of  data  located  at  the  centroid  of  the  2nd  street.  Furthermore  span  wise  mean  velocities  (Vm) 
were  compared  on  the  same  line  of  data.  Samples  of  these  comparisons  for  the  case  of  AOI  =  45°  are 
shown  in  Figure  2.13.  In  this  sample,  calibration  number  7  (purple  line)  meets  these  two  criteria  better  than 
others.  As  discussed  earlier,  the  out-of-plane  velocity  is  measured  in  a  2  mm  thick  layer  which  increases  the 
measurement  sensitivity.  This  issue  is  especially  critical  for  greater  incidence  angles  such  as  AOI  =  45°. 
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Figure  2.9:  Picture  of  the  axisymmetric  jet  and  SPIV  cameras. 


Figure  2.10:  Calibration  target  as  seen  from  both  cameras. 
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Figure  2.11:  Comparison  of  stream  wise  velocity  profiles  for  both  measurement  configurations. 


Figure  2.12:  Deviation  in  %  between  the  mean  streamwise  velocity  measured  in  in-plane  and  out-of-plane 
configurations. 
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Z  /  Height  (mm) 


Vm  vs  Height  rms  velocity  in  Y  direction 


(a) 


(b)  Vrms 


Figure  2.13:  Seven  calibrations,  for  AOI=45°  case. 
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Chapter  3 

Numerical  method 


In  this  chapter,  the  basics  of  the  spectral  element  method  are  presented.  The  spectral  element  discretization 
is  discussed  in  the  context  of  the  Poisson  and  incompressible  Navier-Stokes  equations. 


3.1  Spectral  element  method 

The  Spectral  element  method,  introduced  in  a  1984  paper  by  Patera  [1984],  combines  the  high-order  accu¬ 
racy  of  spectral  methods  with  the  geometric  flexibility  of  the  finite  element  methods.  The  simulation  domain 
(Cl)  is  subdivided  into  spectral  elements  (Clk)  and  the  solution  within  each  spectral  element  is  expressed 
as  a  linear  combination  of  orthogonal  polynomial  basis  functions.  We  define  V^(Clk)  as  the  space  of  all 
polynomials  of  degree  <  TV  in  each  direction  on  domain  Clk.  We  will  restrict  ourselves  to  quadrilateral  and 
hexahedral  spectral  elements. 

A  Vn  approximation  of  the  solution  u(£  1,  £2,  £3)  within  the  standard  hexahedral  region  is  represented 
as  follows: 

N  N  N 

«(£l,£2,6)  =  ^^^UijknN4L(€l)*NjL(&)TTN%L(&)-  C3-1) 

k=0 j= 0 i= 0 

Here,  £1,  £2  and  £3  are  the  local  Cartesian  coordinates,  7r^LL(£)  are  the  one-dimensional  Lagrangian  poly¬ 
nomials  of  order  N  based  on  Gauss-Lobatto-Legendre  (GLL)  nodes  and  utJk  are  the  basis  coefficients.  In 
this  case,  the  total  number  of  basis  functions  is  =  (N  +  l)3. 

3.2  Spectral  element  discretization  of  the  Poisson  equation 

We  first  describe  the  spectral  element  discretization  in  the  context  of  the  Poisson  equation.  The  Poisson 
equation  in  the  domain  Cl  with  homogeneous  Dirichlet  boundary  conditions  is  given  by 

C(u)  =  A 2u  -  f  =  0  in  0  (3.2) 

u  =  0  on  SCI.  (3.3) 

A  typical  spectral  element  mesh  consisting  of  four  spectral  elements  is  shown  in  figure  3.1.  The  spectral 
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Figure  3.1:  A  typical  spectral  element  mesh  consisting  of  4  subdomains. 


element  discretization  is  characterized  by  S  =  ( K ,  N)  where  K  is  the  number  of  spectral  elements  and  N  the 
degree  of  polynomial  approximation  within  each  spectral  element.  We  choose  a  function  space  X$  as 

Xs  =  {<f>\n*  ePN(ttfe)}n?4  (3.4) 

where  Hq  (ft)  is  is  the  space  of  all  functions  which  are  zero  on  the  boundary  and  whose  derivatives  are  square 
integrable  over  the  domain  ft. 

We  assume  that  the  solution  and  the  forcing  function  can  be  accurately  represented  by  a  linear  combina¬ 
tion  of  basis  functions 


Nb 

us(x)  = 

(3.5) 

i= 0 

Nb 

fix)  = 

i= 0 

(3.6) 

where  are  the  basis  (or  trial)  functions  and  Ui,  fi  are  the  expansion  coefficients.  Substituting  the 

approximation  (3.5)  into  equation  (3.2)  produces  the  residual  R  such  that 

C(us)  =  R(us).  (3.7) 

We  introduce  the  Legendre  inner  product  (f,g)  over  the  domain  ft  as 

( f,g)=  [  f{x)g{x)dx.  (3.8) 

The  Method  of  Weighted  Residuals  (MWR)  computes  the  unknown  coefficients  in  equation  (3.5)  by  restrict¬ 
ing  the  residual  to  be  orthogonal  to  a  set  of  test  functions  Vj  (x)  i.e., 

(vj(x),R)=  0.  (3.9) 
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For  our  case  of  the  Galerkin  method,  we  choose  the  trial  functions  to  be  the  same  as  test  functions.  From 
equations  (3.2),  (3.5)  and  (3.7)  we  have  the  variational  form  of  the  Poisson  equation 


Defining  the  stiffness  matrix  K  as 


and  the  mass  matrix  M 


—  /  V'U.Vvdx  =  /  fvdx.. 

Jn  Jn 


■  /  V0i(x)V0j(x)dx 

Jn 

=  /  <Mx)0j(x)dx. 
Jn 


M 


the  expansion  coefficients  for  the  solution  can  be  obtained  by  solving  the  linear  system 


(3.10) 


(3.11) 

(3.12) 


Ku  =  Mf  (3.13) 

3.3  Spectral  element  discretization  of  the  Navier-Stokes  equations 

In  this  section,  we  present  the  spectral  element  discretization  of  the  unsteady  incompressible  Navier-Stokes 
equations.  The  function  spaces  chosen  for  velocity  and  pressure  basis  functions  are  given.  We  then  describe 
the  spatial  discretization  and  temporal  integration  algorithm  for  the  unsteady  Navier-Stokes  equations.  The 
velocity  pressure  decoupling  algorithm  is  discussed.  Finally,  the  multilevel  algorithm  used  for  computation 
of  pressure  is  described. 

3.3.1  Navier-Stokes  equations 

The  Navier-Stokes  equations  governing  incompressible  fluid  flow  in  a  domain  £2  with  boundary  SQ  are  given 
by 


Re( 


dtUi  +  UjdjUi  = 

+ 

■'■'s 

M 

on 

fix  [0  T] , , 

(3.14a) 

diUi  = 

=  0 

on 

fix  [0  T] , , 

(3.14b) 

Ui(t  =  0)  = 

=  u°i 

on 

fi,, 

(3.14c) 

Ui  - 

=  9? 

on 

(5fix>, , 

(3.14d) 

Ui  +  d,UjJj  rij  = 

=  0 

on 

5Q,o, , 

(3.14e) 

Ui{xk)  = 

=  Ui(: 4)  for  xk,  x'k 

on 

50v , 

(3.14f) 

where  Ui(xj),  p(xj )  and  fi(xj)  are  velocity,  pressure  and  body  force  at  each  point  Xj  in  domain  £2.  The 
Reynolds  number  is  Re  =  UL/v ,  where  U  and  L  are  the  characteristic  velocity  and  length  scales  and  v 
is  the  kinematic  viscosity.  Equations  (3.14a)  and  (3.14b)  present  the  momentum  and  continuity  equations 
respectively.  Equations  (3.14c),  (3.14d),  (3.14e)  and  (3.14f)  present  the  initial  condition,  Dirichlet  boundary 
condtion,  outflow  boundary  conditions  and  periodic  boundary  conditions  respectively.  5Qx>,  ££2e>  and  SQp 
represent  the  Dirichlet,  outflow  and  periodic  boundaries  respectively.  In  equation  (3.14f),  x'k  represents  a 
shadow  point  corresponding  to  the  point  Xk  on  the  periodic  boundary. 
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3.3.2  Function  spaces 

Before  proceeding  with  the  spatial  discretization,  it  is  important  to  choose  consistent  function  spaces  for 
velocity  (Xn)  and  pressure  (Yn)  so  that  the  resulting  system  is  stable.  Here,  we  use  the  Vn  —  ^n-2 
formulation  for  velocity  and  pressure.  This  satisfies  the  inf-sup  (Babuska-Brezzi-Ladyzenskaya)  condition 
for  stability.  Tensor  products  of  one  dimensional  Lagrangian  polynomials  based  on  GLL  nodes, 

N  N  N  NVb 

«i(£l,6,&)  =EEE  l)7r£“(&)7r£“(&)  =  £  UiMt  1,6,6)-  (3.15) 

1=0  m=0  n=0  3=0 

are  used  as  basis  functions  for  velocity  within  each  spectral  element  and  C°  continuity  is  enforced  between 
adjacent  spectral  elements.  Here,  NVb  =  (AT  +  l)3  is  the  total  number  of  velocity  basis  functions.  This 
has  certain  important  advantages.  Since  the  interior  basis  functions  are  zero  on  the  boundary  GLL  nodes, 
enforcement  of  Dirichlet  boundary  conditions  and  inter-element  continuity  only  involves  boundary  nodes. 
This  simplifies  the  implementation  of  Dirichlet  boundary  conditions  and  offers  significant  reduction  in  com¬ 
munication  costs.  We  use  ®i(xj)  to  denote  the  global  velocity  basis  functions. 

The  Pn-2  function  space  is  chosen  for  pressure.  Since  inter-element  continuity  is  not  explicitly  en¬ 
forced  for  pressure  between  adjacent  spectral  elements,  a  modal  basis  consisting  of  tensor  products  of  one¬ 
dimensional  Lagrangian  polynomials 


N-2N-2N-2  NPh 

P(€  1,6,6)  =  EEE  PlmnLN,l{£l)LN,m{&)LN,n(t;3)  =  E^ite  !>&,&),  (3.16) 

1=0  rri=0  n= 0  j= 0 

are  used  as  basis  functions  within  each  spectral  element.  Here,  NP^  =  (N  —  1) 3  is  the  total  number  of 
pressure  basis  functions.  We  use  ^i(xj)  to  denote  the  global  pressure  basis  functions. 

3.3.3  Galerkin  projection 

We  now  proceed  to  derive  the  weak  form  of  the  Navier-Stokes  equation  using  the  Galerkin  projection 
method.  The  velocity,  pressure  and  body  force  are  represented  by  a  linear  combination  of  their  global 
basis  functions  as 


NGVb 

iiX'k)  =  ^  ^  & i 


NGPb 


NGVb 


J$j(xk),p(xk)  =  and /*(**)=  E  (3-17) 


3= 0 


3=0 


3=0 


The  Galerkin  method  computes  the  unknown  coefficients  in  (3.14)  by  restricting  the  residual  to  be  orthogonal 
to  the  function  space  used  for  approximating  the  solution.  Thus  the  test  functions  are  the  same  as  the  basis 
functions.  Following  the  variational  procedure  outlined  in  the  earlier  section,  we  obtain 


f  ^k^l  7,  dV  +  f  Uj^rn  'U'ijdV  —  (I)  *&k(TijTlj)dS  f  dj^j^TijdV 

Jn  Jn  J  8Q  Jn 

+  [  Mf\idV, 

Jn 

[  ^kdi$iUi,idV  =  0. 

Jn 


(3.18a) 

(3.18b) 
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We  use  the  convective  form  for  discretizing  the  nonlinear  term.  The  Dirichlet  boundary  conditions  are 
enforced  using  the  lifting  technique,  Ui  =  u]nt  +  uf.  We  can  see  that  the  first  term  on  the  right-hand  side  is 
zero  on  Dirichlet  and  outflow  boundaries  and  can  be  dropped.  These  equations  can  be  rewritten  as 

M~c^t  +  =  ~  Jie^  +  D^' 

DiTLi  —  0, 

where 

M=  f  (j)k(j)kdV,C  =  f  (j)m(j)kdj(j)idV,A  =  f  dkfadkfadV,  and  A  =  [  'ipkdifrdV.  (3.20) 

Here  M,  C  and  A  are  the  mass  matrix,  convective  matrix  and  Laplacian  matrix  respectively.  The  matrix  Di 
contains  the  projection  of  the  derivatives  of  the  velocity  basis  functions  in  the  zth  direction  on  the  pressure 
basis  function. 


(3.19a) 

(3.19b) 


3.3.4  Temporal  discretization 

We  use  a  semi-explicit  scheme  for  temporal  discretization.  The  Stokes  operator,  which  places  a  severe 
restriction  on  time  step  size  if  treated  explicitly,  is  treated  implicitly  using  a  backward  difference  (BDF) 
scheme.  The  nonlinear  convective  term  is  treated  explicitly  using  the  extrapolation  (EX)  scheme.  Discretiz¬ 
ing  equation  (3.19)  using  a  second-order  backward  difference  (BDF2)  scheme  for  the  Stokes  operator  and 
second-order  extrapolation  scheme  (EX2)  we  get 

M  [a0u^+2  -  aiu^+1  -  a2Ui  ]  =  Mf?+2  -  4~Mi+2  +  Dj pn+2 

Re 

-  [(3ou]+1Cu™+1  +  fa u]Cu?]  ,  (3.21a) 

D,u,  =  0,  (3.21b) 


where, 


Using  the  notation 


ao  -  2^’ ai  -  ie tt2  _2 i'P°~ 2’ fji  ~  _L 


u  = 

U2 

,r  = 

h 

h 

,D  = 

'  Di  ' 

d2 

_  ^3  _ 

_  ^3  _ 

.  Dz  . 

M2 


M3 


and  H  =  —  A  +  ao  M. 
Re 


equation  (3.21)  can  be  rewritten  as 


H 

—D 


’  un+2 

r 

_  0  _ 

(3.21c) 


(3.22) 


Here  A,  M  and  H  are  the  global  Laplacian,  mass  and  Helmholtz  matrices  respectively. 
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3.3.5  Pressure  velocity  decoupling 

In  principle,  the  system  of  equations  (3.21)  can  solved  directly  for  pressure  using  the  Uzawa  algorithm.  We 
have 

Sp  =  -DH-1?,  (3.23) 

where 

S  =  DR-xbT .  (3.24) 

For  steady  problems,  the  S  matrix  is  well  conditioned  and  can  be  solved  very  efficiently  using  an  iterative 
method  like  the  preconditioned  conjugate  gradient  method  when  preconditioned  by  the  mass  matrix.  In  this 
case,  while  the  number  of  outer  solves  are  small,  the  H  system  (in  this  case  the  A  system)  needs  to  be  solved 
in  each  iteration. 

However,  for  the  unsteady  Navier-Stokes  equations,  while  the  H  system  is  well  conditioned  and  can 
be  solved  quickly  using  the  diagonally  preconditioned  conjugate  gradient  method,  the  S  matrix  is  very  ill- 
conditioned.  This  significantly  increases  the  number  of  outer  iterations.  This  necessitates  the  use  of  a  more 
efficient  decoupling  algorithm  for  solving  problems  of  engineering  interest. 

In  this  code,  we  use  a  numerical  fractional  step  method  for  decoupling  pressure  and  velocity.  This  is 
discussed  in  detail  in  Couzy  [1995],  Perot  [1993]  and  Fischer  [1997].  The  basic  idea  is  to  compute  the 
pressure  update  only  instead  of  the  full  pressure.  The  system  of  equations  (3.21)  can  be  rewritten  as 


H  -HQDt ' 

un+ 2 

r  +  DT  pn+1  +  res 

—D  0 

pn+ 2  _  pn+ 1 

0 

where  the  splitting  error  is  given  by 

res  =  (I  -  HQ)DT(pn+ 2  -  pn+1).  (3.26) 

Following  Couzy  [1995]  and  Fischer  [1997],  we  choose  Q  =  ^M-1  which  yields  a  splitting  error 

res  =  A—(AM^)Dt (pn+2  -pn+1).  (3.27) 

aoRe 

which  is  second  order  in  time.  In  general,  this  fractional  step  method  can  be  tuned  for  higher  order  of 
accuracy.  Using  Gaussian  elimination,  equation  (3.25)  can  be  reduced  to 

H  — -HM~1Dt 

OL  0 

0  E 

where 

E  =  —DM~1Dt, 
g  =  -DH-1^  +  bTpn+1) 


(3.28a) 

(3.28b) 


r,n+ 2 


-n+ 2 


r  +  DTpn+1  +  res 
9 


The  system  E  is  called  the  consistent  Poisson  operator. 
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3.3.6  Solution  strategy  for  pressure 

The  solution  of  the  Navier-Stokes  equations  requires  two  system  solves  at  each  time  step.  The  first  one 
is  the  solution  of  the  Helmholtz  system.  At  high  Reynolds  number,  the  Helmholtz  operator  is  diagonally 
dominant  and  can  be  solved  in  a  computationally  efficient  way  using  a  Jacobi-preconditioned  conjugate 
gradient  algorithm.  The  most  expensive  step  is  the  solution  of  consistent  Poisson  equation  for  pressure 
given  by 

Ep  =  g.  (3.29) 

We  use  a  two  level  solution  strategy  Fischer  [1996]  for  solving  the  consistent  Poisson  equation.  We  decom¬ 
pose  the  pressure  into  fine  (pf)  and  coarse  (pc)  components.  Figure  3.2  illustrates  a  typical  decomposition 
of  pressure  modes  into  fine  and  coarse  components.  Slowly  varying  pressure  modes  are  chosen  as  coarse 
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Figure  3.2:  Schematic  illustrating  the  decomposition  of  pressure  into  fine  ( pf )  and  coarse  (pc )  modes  for  a 
quadrilateral  spectral  element.  In  this  case,  an  eighth-order  approximation  is  used  for  the  fine  pressure,  and 
a  second-degree  approximation  is  used  for  the  coarse  pressure. 

modes.  The  degree  of  coarse  space  can  be  chosen  at  runtime.  Using  the  decomposition  p  =  pc  +  Pf, 
equation  (3.29)  can  be  rewritten  as 


Eff 

.  Ecf 


’  Pf  ' 

'  9f  ' 

.  P°  . 

.  9c  . 

The  fine  and  coarse  systems  for  pressure  are  given  by 


EfPf  =9f~  EfcEjdc 
EqcPc  =  9c  EcfPf  , 
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(3.30a) 

(3.30b) 


respectively  where  Ef  =  Eff  —  EfcEC(}Ecf.  This  system  is  first  solved  for  pf  and  then  forpc. 


3.3.7  Solution  of  the  fine  pressure  system 


The  fine  pressure  system  is  solved  using  a  block-Jacobi  preconditioned  conjugate  gradient  method.  This 
is  described  in  the  work  of  Couzy  [1995].  We  first  assume  that  the  spectral  element  is  a  rectangle  or  a 
rectangular  parallelpiped  in  three  dimensions.  The  local  consistent  Poisson  operator  Ek  can  be  written  as 


where 


Ek  =  DXM~1DX  +  DyM~lDy  +  DZM~XDTZ 


(3.31) 


Dx  =  <E>  1psWs(f)s  ®  IprWA'r]  , 

Dy  =  [iptWt&t  ®  IpsWsfis  ®  1prWr(f>r}  , 

Dz  =  [iptWt&t  ®  ®  1prWr(j)r\  , 

M_1  =  rrr  I'5*"1  ®  ®  <‘]  - 


(3.32a) 

(3.32b) 

(3.32c) 

(3.32d) 


Here  ^  is  the  matrix  containing  the  one-dimensional  pressure  basis  functions,  </>  is  the  matrix  containing 
the  one-dimensional  velocity  basis  functions  and  <p'  is  the  matrix  containing  the  derivatives  of  one  dimen¬ 
sional  velocity  basis  functions.  These  are  all  evaluated  on  the  Gauss-Legendre  mesh,  w  is  the  diagonal 
matrix  containing  Gauss-Legendre  weights  and  w  is  the  diagonal  matrix  containing  the  Gauss-Lobatto- 
Legendre  weights. 


where 


Ek 


2L  1 


J, 


iEr  +  ^Jt 

Zty 


Jr  +  ‘tE‘ 


Jo  (g)  Jr 


Jr  =  Js  =  Jt  =  lpW(j)W  1 (j)T  wT , 

Er  =  Es  =  Et  =  'i/jwft  u)-1  4>,t wt \pT . 

In  this  form,  Ek  can  be  inverted  using  a  fast  diagonalization  method  (FDM). 


(3.33) 

(3.34a) 

(3.34b) 


4_1  =  (Sz®Sy®Sx)(l^Iz®Iy®Ax  +  lEIz®Ay®Ix  +  lAAz®Iy®Ix)-\sJ®ST®ST)  (3.35) 

Ej  b  y  Zljj; 

where  Sx,  Sy,  Sz,  Ax,  Ay  and  Az  are  the  solutions  of  generalized  eigenvalue  problems 

EXSX  =  JXSXAX ,  EySy  =  JySyAy  and  E ZS z  =  J ZS z Az.  (3.36) 

The  matrix  in  the  middle  is  diagonal  and  can  be  trivially  inverted.  Hence,  the  block  Jacobi  preconditioner  is 
given  by 


rEy 


EpreC 


E1 


E 


K  -I 


(3.37) 
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Figure  3.3:  Schematic  of  a  ’’deformed”  element  (left)  and  the  corresponding 
’’undeformed”  element  (right)  used  for  constructing  the  block  preconditioner  for 
pressure 


In  the  general  case,  the  individual  spectral  elements  are  ’’deformed”  and  the  local  Ek  matrix  cannot  be 
inverted  using  FDM.  In  those  cases,  we  first  construct  a  corresponding  ’’undeformed”  element  with  average 
dimensions  of  the  ’’deformed”  element.  Figure  3.3  shows  a  typical  ’’deformed”  element  and  the  correspond¬ 
ing  ’’undeformed”.  For  example,  in  a  deformed  hexahedral  spectral  element  the  average  seperation  between 
faces  can  be  computed  using 


i:  = 


r  — 

Lt  — 


^2m,n  PmPn[(x Nmn  XQ ran)  ( VNmn  VOmn )  (ZNmn  Z0 mn)  \ 

PmPn 

E l,n  PlPn[(xlNn  ~  XWn Y  +  (VlNn  ~  VlOn)2  +  (ZlNn  ~  Zl0n)2] 

S l,n  PlP™ 

PlPm[(ximN  ~  Xlm())2  +  (VlmN  Vlmo)2  +  (ZtmN  Zlmo)2] 

Sl m  PlPm 


0.5 


(3.38a) 

(3.38b) 

(3.38c) 


Here,  Zr,  ls  and  lt  are  the  average  lengths  in  principal  directions,  N  is  the  order  of  polynomial  approximation 
in  each  direction,  p  are  the  one  dimensional  GLL  weights  and  x,  y  and  z  are  the  cartesian  coordinates  of 
mesh  points. 


3.3.8  Solution  of  the  coarse  pressure  system 

Efficient  solution  of  the  coarse  grid  problem  plays  an  important  part  in  the  scalability  of  the  spectral  element 
solver.  In  this  code,  we  implement  two  coarse  grid  solvers.  The  first  is  a  fast  parallel  direct  solver  described 
in  detail  in  the  work  of  Tufo  [2001].  The  second  is  an  algebraic  multigrid  (AMG)  solver  described  in  the 
work  of  Lottes  [2011].  We  make  two  important  improvements.  The  parallel  direct  solver  is  modified  to 
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support  partition  counts  that  are  not  an  integral  power  of  two.  The  code  for  setting  up  the  AMG  solver  is 
multi-threaded  which  helps  in  reducing  setup  time.  The  parallel  direct  solver  is  used  for  mesh  sizes  less  than 
105  and  the  AMG  solver  is  used  for  larger  mesh  sizes. 

3.3.9  Outflow  boundary  condition  for  turbulent  flows 

For  high  Reynolds  number  flows,  energy  influx  into  the  domain  caused  by  strong  vortices  exiting  at  the 
outflow  boundary  can  create  numerical  instability.  If  we  look  at  the  equation  governing  the  evolution  of 
kinetic  energy, 

TT,  [  \{ukuk)dV  =  v  f  (djUkdjUk)dV  +  [  (, fou^dV 

ot  Jq  z  jn  jn 

+  I  {rijTij  -  ]-ukukrii)uidS  +  f  ( njTij  -  ]-ukukrii)uidS ,  (3.39) 

J dQx>  ^  J  dQo  ^ 

we  can  see  that  the  last  term  on  the  right  hand  side  describing  the  surface  integral  over  dfl  o  can  cause 
numerical  instability  if  rijUj  <  0  anywhere  on  dflo- 

We  implement  the  outflow  boundary  condition  described  in  the  work  of  Dong  et  al.  [2014].  The  idea  is 
to  impose  the  boundary  condition 


rijTij  =  -( ukuk)S0{nkuk)ni  on  Q0,  (3.40) 

so  that  energy  influx  through  the  outflow  boundary  will  not  create  numerical  instability.  Here 

So(nkuk)  =  1  ^1  -  tanh  >  (3-41) 

where  Uo  is  the  characteristic  velocity  scale  and  S  is  a  small  non-dimensional  positive  constant. 

3.3.10  Filter-based  stabilization 

At  high  Reynolds  numbers  a  stabilization  method  is  generally  needed  in  absence  of  a  LES  model.  Here, 
we  use  the  filter-based  stabilization  method  developed  by  Fischer  and  Mullen  [2001].  The  basic  idea  is  to 
interpolate  the  flow  field  onto  a  coarser  mesh  and  interpolate  the  data  back  from  the  coarse  mesh  back  to  the 
fine  mesh.  For  velocity  u  in  Vn,  the  filtered  velocity  ua  is  given  by 

ua  =  (1  —  a)u  +  ijv-i'u,  (3.42) 

where  a  is  the  filter  coefficient.  Here,  In-iu  is  generated  by  interpolating  u  onto  GLL  points  for  Vn-i  and 
interpolating  it  back  onto  GLL  points  for  Vn  •  In  general,  multi-level  filters  can  be  used.  This  procedure  has 
the  advantage  that  inter-element  continuity  is  preserved  and  so  is  spectral  accuracy.  In  general,  the  filter  is 
applied  to  the  flow  field  after  each  time  step. 


3.4  Implementation  details 

In  this  section,  we  present  some  implementation  details  related  to  the  spectral  element  solver  Spec  solve.  The 
code  is  developed  in  C++  and  uses  the  Message  Passing  Interface  (MPI)  for  communication.  It  exploits  the 
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object-oriented  features  of  C++  and  uses  dynamic  memory  allocation  for  optimal  memory  usage.  The  code 
is  relatively  self-contained  and  the  only  external  dependencies  are  the  LAPACK  and  BLAS  libraries.  We 
present  below  a  brief  description  of  some  important  components  of  the  code. 

The  CUBIT  mesh  generator  is  used  for  generating  the  spectral  element  meshes.  Currently  only  quadri¬ 
lateral  and  hexahedral  meshes  are  supported. 

The  program  mesh  .  c  contains  components  for  reading  and  processing  the  spectral  element  mesh.  The 
function  read_mesh  ( )  reads  the  spectral  element  mesh  generated  by  CUBIT.  The  mesh  is  then  scanned 
for  presence  of  periodic  zones.  If  periodic  zones  are  present,  each  face  in  a  periodic  zone  is  matched  with 
appropriate  face  in  the  corresponding  shadow  zone.  It  also  makes  sure  that  the  periodic  and  shadow  faces 
are  properly  aligned.  It  then  matches  periodic  nodes  and  edges  with  corresponding  shadow  nodes  and  edges. 
It  then  assigns  nodes,  edges  and  faces  to  corresponding  cells.  The  data  structures  used  for  nodes,  faces, 
edges,  cells  and  zones  are  defined  in  geom.h.  It  then  proceeds  to  assign  global  id  ( gid )  to  each  degree  of 
freedom.  Figure  3.4  shows  the  spectral  element  mesh  for  a  typical  quadrilateral  spectral  element.  It  numbers 
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Figure  3.4:  Numbering  node  (□),  face  (o)  and  interior  (*)  degrees  of  freedom 
for  a  quadrilateral  spectral  element. 

nodes,  edges  (in  three  dimensions)  and  faces  in  that  order.  The  interior  cell  points  are  numbered  in  the  end 
after  the  entire  skeleton  is  numbered.  For  edges,  faces  and  cells  only  the  gid  range  ([gidfirst,  gidiast  +  1]) 
needs  to  be  stored.  A  weighted  graph  is  then  built  based  on  node  connectivity  information.  The  graph  is 
then  partitioned  using  our  parallel  recursive  spectral  bisection  algorithm.  The  graph  partitioning  algorithm 
is  written  in  graph  .  c.  Figure  3.5  shows  partitioning  of  a  typical  spectral  element  mesh  into  5  partitions. 
This  partitioner  does  not  require  the  number  of  partitions  to  be  an  integral  power  of  two.  This  partition 
information  is  also  used  during  the  setup  stage  by  the  parallel  direct  solver,  which  is  used  to  solve  the  coarse 
grid  problem.  Master  processors  are  assigned  to  each  node,  edge,  face  and  cell.  It  then  assigns  boundary 
conditions  to  each  degree  of  freedom. 

The  program  solver  .  c  initializes  the  ’’Solver”  object.  The  Solver  class  is  defined  in  solver.h.  This  pro¬ 
gram  distributes  mesh  data  to  each  processor  based  on  partition  information.  In  addition,  mesh  information 
related  to  cells  on  processor  boundaries  are  distributed  redundantly  to  all  processors  involved. 

The  program  bas  i  s  .  c  generates  the  one-dimensional  basis  functions,  quadrature  points  and  quadrature 
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Figure  3.5:  A  two-dimensional  spectral  element  mesh  partitioned  into  5  parts  using 
recursive  spectral  bisection. 


weights.  Figure  3.6  shows  the  one  dimensional  velocity  and  pressure  basis  functions  associated  with  a 
quadrilateral  spectral  element.  Figure  3.7  shows  three  different  meshes  for  evaluating  various  integrals. 
The  pressure  mesh  is  used  for  evaluating  integrals  involving  pressure  basis  functions.  Mesh  2  is  used  for 
de-aliased  integration  of  nonlinear  terms.  Mesh  1  is  used  for  computing  all  other  integrals.  Mesh  1  and 
2  used  GLL  (Gauss-Lobatto-Legendre)  quadrature  of  appropriate  order  whereas  the  pressure  mesh  uses 
LG(Legendre-Gauss)  quadrature. 

The  program  matrix,  c  defines  various  matrix  classes  and  related  functions  used  by  the  program. 
Serial  dense  matrices  and  parallel  sparse  matrices  are  supported.  A  sparse  matrix-vector  product  is  custom 
coded  whereas  LAPACK  and  BLAS  are  used  for  dense  matrix  operations.  The  file  lpk_blas  .  c  provides 
the  interface  to  external  LAPACK  and  BLAS  libraries. 

The  program  map2stdel .  c  contains  code  for  generating  the  internal  spectral  element  mesh.  First, 
nodes  on  the  exterior  of  each  cell  are  mapped  to  appropriate  locations  and  a  smooth  interior  nodal  distribu¬ 
tion  is  obtained  by  solving  a  Poisson  equation  with  exterior  node  distribution  as  boundary  condition.  It  also 
computes  Jacobian’s  and  other  entities  needed  for  local  elemental  operations  like  integration  and  differenti¬ 
ation. 

The  program  dss  .  c  contains  functions  for  setting  and  effecting  direct  stiffness  summation.  The  pro¬ 
gram  tens  .  c  contains  functions  that  evaluate  the  effect  of  an  operator  on  a  vector.  In  this  code,  none  of  the 
operators  are  explicitly  constructed.  For  quadrilateral  and  hexahedral  spectral  elements,  these  operators  can 
be  written  as  a  tensor  product  (see  Deville  et  al.  [2002]).  For  a  d-dimensional  problem,  the  sum  factorization 
technique  can  be  used  to  evaluate  these  tensor  product- vector  multiplications  in  0{nd+1)  operations  instead 
of  the  naive  0(n2d)  operation  count  for  a  matrix  of  order  n  . 

The  GMRES  (generalized  minimal  residual)  solver  for  pressure  and  the  PCG  (preconditioned  conjugate 
gradient)  solver  for  the  Helmholtz  system  are  written  in  gmres  .  c  and  peg .  c  respectively.  The  fractional 
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(a)  Velocity  basis 

Figure  3.6:  One-dimensional  basis  functions  for  T$ 


(b)  Pressure  basis 

—  V3  spectral  element  formulation. 


step  method  is  implemented  in  navier_stokes.c.  The  program  coarse .  c  contains  functions  for 
building  the  coarse  pressure  matrix.  The  direct  parallel  solver  for  the  coarse  pressure  problem  is  coded  in 
xx t .  c.  The  algebraic  multigrid  solver  is  coded  in  amg .  c. 

For  meshes  containing  more  than  100000  spectral  elements,  AMG  solver  needs  to  be  used.  A  multi¬ 
threaded  solver  setup_amg  is  built  for  generating  data  needed  for  the  AMG  solver.  In  this  case,  first  the 
coarse  pressure  matrix  is  built  and  written  to  a  file.  setup_amg  is  called  to  set  up  data  needed  for  the  AMG 
solver  before  running  the  spectral  element  solver. 

Finally,  the  program  export .  c  contains  functions  that  perform  input-output  (IO)  operations.  We  use 
MPI-IO  for  optimal  IO  performance.  Visit  and  Tecplot  are  used  for  visualization  of  results. 


3.5  Algebraic  multigrid  solver 

In  this  section,  we  describe  in  detail  the  algebraic  multigrid  solver  used  for  the  solution  of  the  coarse  problem. 

3.5.1  Introduction 

We  consider  the  solution  of  the  linear  system 


Ax  =  b,  (3.43) 

where  A  is  a  sparse  symmetric  positive  definite  n  x  n  matrix  using  fixed  point  iteration  method.  Using  the 

initial  guess  xinit  and  a  symmetric  positive  definite  smoother  C  =  B-1,  the  iterative  solution  method  can 
be  written  as 

x°  =  xinit,  (3.44) 

xm+l  =  xrn  +  C(b_Axm^  (3.45) 
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(a)  Mesh  1  (GLL  (N)) 


(b)  Mesh  2  (GLL  (3N/2))  (c)  Pressure  mesh  (LG  (N-2)) 


Figure  3.7:  Quadrature  points  for  velocity  mesh. 


with  the  corresponding  error 


em+i  =  (/  _  CA)em  =  (I  -  CA)me°.  (3.46) 

In  general  C  can  be  allowed  to  be  singular  as  is  the  case  with  the  coarse-grid  correction  step  of  the  algebraic 
multigrid  method.  The  convergence  of  this  algorithm  is  controlled  by  the  spectral  radius  of  the  matrix 
(/  —  C A).  We  consider  the  generalized  eigenvalue  problem 


Azk  =  \kBzk,  (3.47) 

where  zk  are  the  orthonormal  eigenvectors,  Xk  are  the  eigenvalues  and  B  is  scaled  so  that  0<Ai<A2< 
*  *  •  <  An  <  1.  Using  em  =  ^ k  ekmzk ,  it  follows  from  equation  (3.46)  that 

ekm  =  (1  -  Afc)mefc°,  1  <  k  <  n.  (3.48) 

We  can  see  that  the  smoother  is  effective  at  reducing  the  high- wavenumber  error  components  whereas  it  is 
slow  at  reducing  the  low- wavenumber  error  components. 

The  basic  idea  of  the  multigrid  method  is  to  combine  a  smoother,  which  is  effective  at  reducing  the 
high- wavenumber  error  with  a  coarse-grid  correction  to  target  the  low- wavenumber  error  components.  A 
typical  multigrid  algorithm  consists  of  three  steps.  The  first  step  is  smoothing,  which  is  effective  at  reducing 
the  high-wavenumber  error.  The  second  step  involves  restriction  of  the  residual  to  the  coarse  level  and 
computation  of  the  coarse-level  correction.  The  third  step  involves  interpolation  of  the  coarse  correction  to 
the  fine  level.  This  could  be  accompanied  by  further  post  smoothing.  An  L-level  multigrid  algorithm  can 
be  completely  characterized  by  the  coarsening  algorithm,  smoother  C/,  and  the  prolongation  operator  Pi  at 
each  level.  Defining  the  two-level  multigrid  error  propagation  matrix  at  the  coarsest  level  Mtg ?l,  as 

Etg,L  =  ( Il  ~  CLAL)mu’L{IL  -  PLA^PlAL){IL  -  CLAL)md’L ,  (3.49) 

and  the  two-level  multigrid  error  progagation  matrix  at  level  fc,  k  <  L  ,  as 

Etg,k  =  (Ik  ~  CkAk)m^(Ik  -Pk(I-  E%tk+1)Akhpk  Ak)(h  ~  CkAk)md’\  (3.50) 
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a  multigrid  iteration  can  be  written  as 


x-xi+1  =  £ms(x-Xi),  (3.51) 

where  the  multigrid  error  propagation  matrix  Emg  =  Etg  o,  can  be  computed  recursively.  Here  Aq  =  A, 
Ak+1  =  PfAkPk  is  the  coarse  matrix  at  each  level,  ra^,  rau+  are  the  number  of  pre-  and  post- smoothing 
iterations,  respectively,  at  each  level,  and  7^  is  the  number  of  lower-level  multigrid  iterations  made  at  each 
level.  The  convergence  rate  of  the  multigrid  method  is  determined  by  the  spectral  radius  of  Emg ,  p(Emg). 
While  convergence  analysis  of  the  entire  multigrid  algorithm  is  quite  complicated,  algorithms  built  from 
efficiently  constructed  two-level  multigrid  algoritms  with  small  p(Etg^)  results  in  multigrid  algorithm  with 
good  overall  convergence  properties. 


3.5.2  Two-level  multigrid 


We  now  look  at  the  two-level  multigrid  iteration.  We  start  by  decomposing  the  variables  into  nc  coarse 
variables  (”C-variables”)  and  rtf  fine  variables  (”F- variables”).  The  matrix  A  and  the  prolongation  matrix 
P  can  be  written  in  block  form  as 


Aff 

Afc 

p  — 

'W' 

As 

ACC 

5  1  — 

I 

Am 


Here,  W  is  the  rtf  x  nc  matrix  consisting  of  interpolation  weights. 
A  two-level  multigrid  iteration  can  be  written  as 


(3.52) 


'■2  +  1 


=  -Etg(x-Xi)f 


(3.53) 

(3.54) 


where 

Etg  =  {I-  BA)mu(I  -  PApPTA) 

is  the  two-level  error  propagation  matrix.  Here  B  is  a  smoother  and  Ac  =  PTAP  is  the  coarse  matrix.  In 
this  case,  we  only  use  a  post-smoother.  The  convergence  rate  of  this  two-level  method  is  determined  by  the 
spectral  radius  of  Etg,  p(Etg). 

Following  Lottes  [2011]  using  the  hierarchial  basis 


7  W 

, A = Tt AT = 

Aff 

Afc 

I 

-Acf 

Ac_ 

,  b  =  TTb,  x  =  Tx, 


(3.55) 


the  linear  system  (3.43)  can  be  transformed  to 


Aff 

Afc 

V 

V 

-Acf 

Ac. 

A. 

be. 

(3.56) 


where  Afc  =  AffW  +  Afc  and  Ac  =  PTAP.  The  prolongation  and  smoother  matrices  transform  to 


P 

II 

1 

ba 

II 

B 

=  T~lBT~ 

0 
Ia 

Bjf  Bfc 
Be  f  Bcc_ 

Bff  —  BfcWT  -  WBcf  +  WBccWt  Bfc-  WBc 
Bcf  -  BCCWT  Bcc 


(3.57) 


(3.58) 
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respectively,  where 


T~ 1  = 


I 

0 


-w 

I 


(3.59) 


The  choice  of  weights  W  =  A^Afc  is  ideal  and  decouples  the  fine  and  coarse  problem.  However  in  this 
case  the  weight  matrix  W  is  full  and  is  avoided.  The  error  propagation  matrix  is 


Etg  =  (I  -  BA)mu(I  -  PAc1PtA)  =  T-'EtgT. 


(3.60) 


We  can  see  that  Etg  and  Etg  have  identical  spectra.  The  quality  of  the  interpolation  is  defined  by 

7=||  A1/? FA-1'2  ||,  (3.61) 

where  F  =  W  —  (—Aj^Afc).  The  coarse  system,  Acxc  =  bc,  is  solved  first.  This  needs  to  be  solved 
recursively  if  it  is  not  the  coarsest  level.  The  fine  system,  Affitf  =  bf  —  Afcxc ,  is  solved  next. 


Error  estimates 

'Sf  ol 

Ac\  ’ 

where 

Sf  =  A/7  —  AfcAc  Afc, 


Q  = 


~Ac  lAJc 


, A = QtAQ = 


(3.62) 

(3.63) 


P 

B 


Q-'P 


0 

I  ’ 


q~xbq~t 


Bff  Bfc 
Bcf  Bcc 


.  Bff  „ 

LBff  +  Bcf 


BffLT  +  Bfc 

BBffLT  +  LBfc  +  BcfLT  +  Bcc 


where  L  =  Ac  1  Ajc 


Q~ 


a-^AJc  i. 


Etg  =  (I-  BA)mpI  -  PAc1PtA)  =  Q-'EtgQ. 


I-PA~1PtA  =  I  - 


0 

\Sf 

7 

47 

Ac_ 

0 

Etg  =  (I-  BA) 


7 

1  -  Bjfsf  0 

o_ 

—Bcf  0 

Bff  =  Bff=[l  —W]  B  [I  -W]T 


(3.64) 

(3.65) 

(3.66) 

(3.67) 

(3.68) 

(3.69) 

(3.70) 
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(3.71) 

/  _  II  A  2  p  / 1  I  II  _  ^lln 

/- \\Af  ft  Ac  \\2  -  sup 

V#0  v  Uo 

(3.72) 

p(Etg)  <  1  —  (1  —  72)(1  —  Pf) 

(3.73) 

S 

’ — 1  |CN 

Q 

s 

s 

a 

V 

II 

s 

se 

(3.74) 

7  =  114/^^112 

(3.75) 

Pf  =  Pi1  - 

(3.76) 

-Aj}Afc),Bff=[l  —W]  B  [I  -W]T 

(3.77) 

We  now  describe  the  smoothing,  coarsening  and  interpolation  steps  in  detail. 

Smoother 

We  use  a  smoother  of  the  form 

for  solving  the  fine  system. 

Coarsening 

X  =  I~  D~JAffD~} 

Kf  =  K(I-X) 


B  = 


Bff  0 
0  0 


(3.78) 


n  =  eJ\X\l 

1  H-  V max  t 

Kf  =  — - when  Vmax  <  1 

-t  ^  max 


(3.79) 

(3.80) 


Interpolation 


=  ii4^c|ii2 

\\a)sfd\\f 

minimize  tr(PT AP)  subject  to  Wu  =  —A^Afc u 

(3.81) 

X^RjiRiAffRfy'Ri, 

(3.82) 
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G  2  —  -2Q(  Afc&i  +  UiX), 


(3.83) 


where 


XX  =  -AffAfc u  -  Y^UiXi(-Afcei),X  =  jPufXj. 

i= 1  i=l 

-  ( DAff  +  AffD)  «  X”1,!)  =  diag(d~1),  dj  =  Y,uleJRiRiej- 


i=  1 
x-l\ 


7=  || AffAfcA~%  «  ||D//M/c£>c-1||2, 

(r+p)  , - 

p||2P  <  max  ,w(r)  =  (i?-BT)rl,Cj  =  y  wf } /wf}, 

7  PS  ||D7/raCl2i/cD“1||2  <  Cmax,Cmax  =  max  q. 

JJ  l<i<nc 


(3.84) 

(3.85) 

(3.86) 

(3.87) 

(3.88) 


Table  3.1:  AMG  setup  for  flow  over  2d  circular  cylinder  simulation,  tconv  =  0.1, 
ltarg  =  0.226. 


Level 

N 

Nf 

Nc 

Nc/N 

PS 

m 

7 

nnz(W)/nc 

nnz(Aff)/nf 

1 

199 

125 

74 

0.371 

0.564 

4 

0.201 

7.324 

4.120 

2 

74 

53 

21 

0.283 

0.678 

4 

0.315 

11.142 

11.867 

3 

21 

17 

4 

0.190 

0.656 

4 

0.197 

13.750 

15.941 

4 

4 

3 

1 

0.250 

0.594 

4 

0.000 

3.000 

3.000 

Table  3.2:  AMG  setup  for  flow  over  wall  mounted  cube  simulation,  tconv  =  0.5, 
'"ftarg  =  0.541. 


Level 

N 

Nf 

iVc 

Nc/N 

Pf 

m 

7 

nnz(W)/nc 

nnz(Aff)/nf 

1 

18362 

7972 

10390 

0.565 

0.804 

4 

- 

20.077 

8.865 

2 

10390 

6040 

4350 

0.418 

0.614 

4 

- 

50.152 

144.717 

3 

4350 

3000 

1350 

0.310 

0.677 

4 

- 

119.565 

750.170 

4 

1350 

1085 

265 

0.196 

0.784 

5 

0.292 

227.566 

818.321 

5 

265 

227 

38 

0.143 

0.774 

4 

0.371 

127.210 

227.000 

6 

38 

34 

4 

0.105 

0.776 

4 

0.173 

31.750 

34.000 

7 

4 

3 

1 

0.250 

0.421 

3 

0.000 

3.000 

3.000 
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Table  3.3:  AMG  setup  for  3d  urban  boundary  layer  array  simulation,  tconv  =  0.5, 
7 targ  0.541. 


Level 

N 

Nf 

Nc 

Nc/N 

pf 

m 

7 

nnz(W)/nc 

nnz(Aff)/nf 

1 

164604 

81896 

82708 

0.502 

0.844 

3 

- 

11.856 

9.719 

2 

82708 

48966 

33742 

0.407 

0.707 

3 

- 

39.205 

77.497 

3 

33742 

23990 

9752 

0.289 

0.746 

3 

- 

101.838 

485.489 

4 

9752 

7413 

2339 

0.239 

0.758 

3 

- 

158.084 

1084.056 

5 

2339 

1902 

437 

0.186 

0.777 

3 

0.531 

202.924 

1102.716 

6 

437 

366 

71 

0.162 

0.742 

2 

0.497 

107.112 

366.000 

7 

71 

61 

10 

0.140 

0.711 

2 

0.568 

27.900 

61.000 

8 

10 

9 

1 

0.100 

0.452 

2 

0.000 

9.000 

9.000 
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Chapter  4 


Experimental  results 


4.1  Flow  characteristic  study  tools. 

Two-dimensional  contour  plots  are  used  to  present  and  compare  the  flow  characteristics  such  as  TKE, 
Reynolds  stresses  and  mean  velocities  and  are  used  to  study  the  three-dimensional  behavior  of  the  turbu¬ 
lent  flow  in  our  regions  of  interest.  Figure  4.1  shows  a  sample  for  a  group  of  x  —  z  slices  at  different 
span  wise  positions  with  the  color  indicating  the  mean  streamwise  velocity.  Horizontal  x  —  y  slices  (Sz) 
and  vertical  x  —  z  slices  (Sy)  will  be  used  in  the  discussion  to  depict  various  flow  quantities;  specifically,  a 
combination  of  one  Sy,  and  two  Sz  at  different  heights  will  be  used. 

According  to  Stull  [1988]  the  region  of  the  boundary  layer  directly  above  the  mean  height  of  the  blocks 
is  denoted  as  the  ‘wake  layer’  and  the  region  beneath  is  called  the  ‘urban  canopy’.  The  region  of  large 
initial  turbulence  in  the  wake  layer  of  the  urban  array  is  mostly  due  to  the  flow  separation  occurring  at 
the  leading  edge  of  the  upstream  block  (Martinuzzi  and  Havel  [2000]).  In  this  study  each  street  area  is 
divided  into  two  main  sub-regions;  the  ‘lower  layer’  starts  from  ground  level  up  to  the  mid-height  of  a  block 
(0  <  z/H  <  0.5)  and  the  ‘upper  layer’  includes  the  upper  half  of  each  street  (0.5  <  z/H  <  1).  The  region 
between  two  adjacent  spanwise  blocks  is  called  the  intersection.  Figure  4.2(a)  shows  horizontal  planes  in  the 
lower  layer  (red)  at  z/H  =  0.25  and  upper  layer  (red)  at  z/H  =  0.85  as  well  as  the  intersection  (purple). 
In  addition  a  vertical  slice  (Sy)  is  chosen  at  the  mid-span  plane  of  each  street  (y/W  =  0)  shown  as  the 
green  slice  in  Figure  4.2b.  The  combination  of  these  three  slices  are  used  to  investigate  the  complicated 
phenomena  in  the  urban  boundary  layer  investigated  in  this  study.  Four  flow  quantities  are  studied  in  each 
of  these  three  slices:  mean  velocity  (J7m,  Vm  and  Wm),  turbulent  kinetic  energy  (TKE)  and  its  components 
( u /2,  u/2and  w'2)  and  Reynolds  stress  components  (—u'vf,  —u'w' and  —v'w'). 

4.2  Results 

In  this  section,  different  flow  characteristics  are  studied  for  four  angle  of  incidence  (AOI)  values:  0°,  15°,  30° 
and  45°.  In  order  to  simplify  the  evaluations,  a  street- to- street  comparison  is  run  for  each  angle  of  incidence 
(AOI)  condition,  and  then,  the  effect  of  AOI  on  different  streets  is  studied.  The  effects  are  discussed  both 
in  the  lower  layer  (at  z/H  =  0.25)  and  upper  layer  (z/H  =  0.85)  slices  as  well  as  the  center  vertical  slice 
(y/W  =  0) .  ' 
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Figure  4.1:  Sample  of  a  group  of  vertical  slices  showing  streamwise  mean  velocity. 


(a)  Geometrical  categorization  of  street  region 


Figure  4.2:  Regions  under  study. 
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4.2.1  Mean  velocity 
Contours  in  Sy  slices 

Stream  wise  mean  velocity.  All  mean  velocity  magnitudes  are  normalized  by  Uh  ,  which  is  the  stream  wise 
component  of  the  mean  free  stream  flow  velocity  at  the  height  of  a  building  ( H  =  50  mm).  For  AOI  =  0°, 
the  Sy  slice  at  mid- span  (y/W  =  0)  for  Um  shows  the  high-speed  free- stream  flow  at  the  upper  roof  level 
(z/H  >  1.3)  in  all  streets  (Figure  4.3).  A  strong  velocity  gradient  (^)  is  detected  in  a  thin  layer  near 
the  roof  level  (0.8  <  z/H  <  1.2)  for  all  streets.  At  lower  values  (z/H  <  0.8),  an  approximately  zero 
streamwise  mean  velocity  is  observed  close  to  both  the  upstream  and  downstream  buildings’,  with  a  slightly 
negative  values  observed  in  the  central  part  of  the  streets.  The  contours  shown  in  4.3  show  that  for  AOI  =15° 
the  streamwise  velocity  distribution  on  the  x  —  z  mid- span  plane  for  all  streets  is  very  similar  to  AOI  =  0° 
except  for  the  upper  roof  level  of  the  2nd  ,  3rd  and  4th  streets  in  which  the  shear  layer  is  observed  in  a  larger 
region  compared  to  AOI  =  0°.  The  changes  of  the  distribution  of  the  streamwise  component  of  the  mean 
velocity  is  considerable  in  transition  from  AOI  =  15°  to  AOI  =  30°.  The  Sy  slice  for  Um/Un  shows  that 
for  AOI  =  30°  the  streamwise  velocity  is  always  positive  in  the  entire  street  at  the  mid- span  of  the  streets. 
Moreover,  there  is  a  region  of  slightly  positive  streamwise  velocity  in  the  downstream  half  of  the  1st  and  2nd 
streets.  Also,  the  thickness  of  the  shear  layer  at  the  upper-roof  level  is  thicker  than  AOI  =15°,  and  grows 
in  transition  from  the  1st  to  the  4th  street.  The  thickness  of  the  shear  layer  is  even  larger  for  AOI  =  45°:  it 
is  detectable  in  the  range  of  0.8  <  z/H  <  1.4.  Also,  the  streamwise  velocity  magnitudes  at  the  upper-roof 
level  are  considerably  smaller  for  AOI  =  45°  as  compared  to  the  smaller  AOI  conditions  investigated. 
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(a)  AOI  =  0° 


(c)  AOI  =  30° 


Figure  4.3:  Um/Un  in  the  Sy  plane. 


10.8 
0.6 
0.4 
0.2 
0 

-0.2 


(b)  AOI  =  15° 


(d)  AOI  =  45° 


Span  wise  mean  velocity.  The  mean  span  wise  velocity  distributions  at  the  mid-span  x—z  plane  (y/W  =  0) 
are  shown  in  Figure  4.4.  For  AOI  =  0°,  zero  values  for  Vm/Un  are  expected  based  on  symmetry  and  a 
significant  change  is  observed  with  increasing  AOI.  Higher  magnitudes  of  Vm  are  observed  in  the  streets 
near  the  downstream  building.  The  magnitudes  are  greater  in  the  1st  street  as  compared  to  the  2nd,  3rd 
and  4th.  By  increasing  the  AOI  to  30°  and  45°,  the  spatial  extent  of  the  non-zero  values  of  the  velocity 
increase;  however,  the  magnitudes  show  a  gradual  decrease.  For  the  upper-roof  region  (1.0  <  z/H  <  1.5) 
the  magnitude  of  the  span  wise  velocity  increases  from  AOI  =  0°  to  AOI  =  45°. 
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Figure  4.4:  Vm/Un  in  the  Sy  plane. 


Wall-normal  mean  velocity.  The  Sy  slices  of  wall-normal  mean  velocity,  Wm/Un ,  contours  (Figure 
4.5)  for  AOI  =  0°  show  that  there  are  two  main  regions  of  wall-normal  flow  in  all  four  streets:  a  negative 
(downward)  flow  close  to  the  wall  of  the  downstream  building,  and  a  positive  (upward)  flow  close  to  the 
wall  of  the  upstream  building.  For  all  four  AOI  values,  the  effect  of  the  relatively  strong  flow  separation  off 
the  first  block  is  observed  (see  Martinuzzi  and  Havel  [2004]).  In  transition  from  the  1st  to  the  4th  street, 
the  mean  velocity  magnitude  of  the  downward  flow  near  the  downstream  building  increases  but  the  mean 
velocity  values  of  the  upward  flow  near  the  upstream  building  are  decreasing.  This  combination  of  the 
downward  and  upward  flow  regions  within  the  streets  is  evidence  of  a  single  recirculation  region  within  the 
street  and  is  expected  for  the  skimming  flow  regime. 
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(a)  AOI  =  0° 


(b)  AOI  =  15° 


(d)  AOI  =  45° 


Figure  4.5:  Wm/Un  in  the  Sy  plane. 
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Contours  in  Sz  slices 


Streamwise  mean  velocity.  The  Sz  slices  at  wall-normal  heights  of  z/H  =  0.25  and  z/H  =  0.85  are 
shown  in  Figure  4.6.  For  AOI  =  0°,  symmetric  flow  is  observed  with  higher-speed  Um  values  in  the  intersec¬ 
tions  as  expected.  Um  values  in  the  intersection  are  higher  in  the  1st  street  and  decrease  as  the  flow  moves 
downstream.  There  is  a  region  of  zero  Um  values  close  to  the  upstream  and  downstream  buildings’  walls 
and  was  observed  for  the  Sy  slices  discussed  above.  Also  there  is  an  area  of  reverse  flow  at  the  center  of  the 
streets.  This  reverse  flow  together  with  the  positive  streamwise  flow  on  the  sides  shows  the  existence  of  the 
legs  of  the  classical  arch  vortex  within  the  street  region.  The  effect  of  AOI  is  clearly  seen  in  the  figures  (note 
that  the  direction  of  the  incoming  flow  is  shown  by  the  arrow  on  the  first  block).  This  changing  distribution 
of  mean  flow  results  in  a  shift  of  the  classical  arch  vortex  and  will  be  discussed  in  more  detail  later. 


(a)  AOI  =  0°,  z/H  =  0.25 


(b)  AOI  =  0°,  z/H  =  0.85 


0.8 

0.6 

0.4 

0.2 

0 

-0.2 


(c)  AOI  =  15°,  z/H  =  0.25 


(d)  AOI  =  15°,  z/H  =  0.85 


Figure  4.6:  Streamwise  mean  velocity  normalized  by  Uh  in  Sz  slices  at  z/H  =  0.25  and  z/H  =  0.85. 


Span  wise  and  Wall-Normal  mean  velocity.  The  contours  of  the  normalized  span  wise  (Vm  /Uh)  and  Wall 
normal  (Wm/Un)  components  of  the  mean  velocity  are  shown  in  Figures  4.7  and  4.8  respectively  for  two 
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horizontal  planes  in  all  four  streets.  Fairly  symmetric  results  for  AOI  =  0°  are  observed  along  with  a  strong 
dependence  on  AOI.  For  example  the  spanwise  channelling  of  the  flow  increases  with  increasing  AOL  is 
clearly  seen  in  Figure  4.7.  Also  for  the  wall  normal  flow  (4.8)  increases  in  the  upward  flow  for  y/W  <00 
and  downward  flow  for  y/W  >  0  are  seen  with  increasing  AOI:  this  effect  is  particularly  apparent  higher  in 
the  street  at  z/H  =  0.85. 


(e)  AOI  =  30°,  z/H  =  0.25 


(f)  AOI  =  30°,  z/H  =  0.85 
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(g)  AOI  =  45°,  z/H  =  0.25 
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(h)  AOI  =  45°,  z/H  =  0.85 


Figure  4.7:  Spanwise  mean  velocity  normalized  by  Uh  in  Sz  slices  at  z/H  =  0.25  and  z/H  =  0.85. 
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(a)  AOI  =  0°,  z/H  =  0.25 


(b)  AOI  =  0°,  z/H  =  0.85 


(d)  AOI  =  15°,  z/H  =  0.85 


(e)  AOI  =  30 °,z/H  =  0.25 


(f)  AOI  =  30 °,z/H  =  0.85 


Figure  4.8:  Wall  normal  mean  velocity  normalized  by  Uh  in  Sz  slices  at  z/H  =  0.25  and  z/H  =  0.85 
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4.2.2  Turbulent  kinetic  energy 


Turbulent  kinetic  energy  (TKE),  as  defined  by  Equation  4.1,  will  be  normalized  by  Ujj  and  will  be  denoted 


as  TKE//. 


TKE  =  ^  (u'2  +  v'2  +  w 


(4.1) 


where  uf,  v' ,  and  u/  are  the  fluctuating  components  of  the  velocity  as  defined  using  the  classic  Reynolds 
decomposition  with  the  overbar  representing  time  averaging.  TKE//  is  therefore: 


TKE// 


TKE 

UH 


Figures  4.9  and  4.10  show  the  dependence  of  the  TKE//  with  AOI  in  a  vertical  slice  at  midspan  in 
the  urban  array  at  y/W  =  0  for  the  four  AOI  cases  investigated.  Also  included  in  the  figures  are  the 
contributions  of  each  of  the  three  components  of  the  fluctuating  velocities  u '2,  \  i/2and  \  w'2)  to  the 
overall  TKE  with  color  indicating  the  magnitude  of  the  normalized  values.  It  is  instructive  to  look  at  the 
relative  magnitudes  of  the  component  contributions  in  the  various  regions  of  the  urban  array.  For  example, 
the  high  TKE  levels  just  upstream  of  the  downstream  block  are  mainly  due  to  the  span  wise  component  of 
the  fluctuation  velocity  for  AOI  =  0°  and  15°  as  seen  in  Figure  4.9.  In  the  separating  shear  layer  near  z/H  = 
1  the  high  TKE  values  are  due  to  both  the  streamwise  and  span  wise  components  of  the  fluctuation  velocities. 
The  wall  normal  component  w'2/Ufj)  contributes  very  little  to  the  TKE  in  the  street.  As  AOI  increases  to 
30°  and  45°  (Figure  4.10)  the  overall  levels  of  normalized  TKE  decrease  because  of  the  strong  channelling 
of  the  overall  flow  but  it  is  still  observed  that  the  main  contribution  to  the  shear  layer  at  z/H  =  1  is  due  to  the 
streamwise  and  spanwise  components  and  that  the  wall-normal  component  does  not  contribute  significantly 
to  the  overall  TKE  at  the  midspan  of  the  street. 

The  horizontal  Sz  contour  slices  (at  z/H  =  0.25)  of  TKE//  and  the  relative  contributions  of  the  fluctu¬ 
ating  velocity  components  (|  u'2,  \  u'2and  \  w'2)  are  provided  in  Figures  4.11  and  4.12  for  the  four  AOI 
values  investigated.  A  slight  asymmetry  is  observed  in  the  TKE//  distribution  in  street  1  for  AOI  =  0°.  Every 
effort  was  made  to  place  the  urban  array  in  the  wind  tunnel  at  an  AOI  =  0° ;  however,  a  slight  asymmetry  is 
observed.  Monnier  et  al.  [2010]  also  reported  a  similar  issue:  they  found  that  for  even  a  very  small  deviation 
in  AOI  (  =  0.5°)  asymmetry  in  flow  characteristics  can  be  observed.  The  main  contribution  to  the  TKE  in 
the  intersections  is  from  \  u'2,  for  AOI  =  0°  and  the  contribution  to  TKE  in  the  region  just  upstream  of 
the  downstream  block  is  from  the  spanwise  component  \  v'2.  Again  the  contribution  of  the  wall  normal 
component  \  w'2,  to  the  TKE  is  negligible  in  the  entire  street.  As  the  AOI  is  increased  to  30°  and  45°  the 
normalized  TKE  values  decrease  in  the  street  overall.  The  contribution  to  the  TKE  in  the  intersections  are 
not  just  due  to  -  u'2,  but  also  to  the  spanwise  component  \  v'2. 


4.2.3  Reynolds  Shear  Stress 

The  results  for  normalized  —u'w' and  —  u'v'are  shown  in  Figures  4.13  and  4.14  respectively.  For  the  stream- 
wise  and  wall-normal  stress  we  choose  to  show  vertical  slices  at  mid  span  in  Figure  4.13.  For  all  four  AOI 
values  the  largest  magnitudes  are  due  to  the  separating  shear  layer  from  the  top  surface  of  the  block.  Again 
due  to  the  relative  strong  separation  for  street  1  and  compared  with  the  downstream  streets  the  magnitudes 
of  —u' w' decrease  with  downstream  distance.  The  normalized  —  'uVvalues  are  provided  in  horizontal  planes 
at  z/H  =  0.25  in  Figure  4.14.  The  expected  distribution  at  AOI  =  0°  is  observed;  namely  the  sign  of  —  u'v'is 
positive  in  the  region  0  <  y/H  <  1  and  negative  for  -1  <  y/H  <  0.  As  AOI  increases  the  stream wise- 
spanwise  correlation  distribution  is  significantly  affected  by  the  channeling  effect  due  to  the  non  zero  AOI 
values. 


48 


(a)  AOI  =  0°,  TKEh  (b)  AOI  =  15°,  TKEW 


(g)  AOI  =  0°,  \  w'2/Ujj 


(h)  AOI  =  15°,  \  wl2/Ujj 


Figure  4.9:  TKE#  and  the  contribution  of  fluctuating  velocity  components  normalized  by  Ufr ,  in  Sy  slices 
at  y/W  =  0  for  AOI  =  0°  and  AOI  =  15°. 
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(a)  AOI  =  30°,  TKEH  (b)  AOI  =  45°,  TKEW 


1.5 


(g)  AOI  =  30°,  |  w'Z/Ufj 


Figure  4.10:  TKE//  and  the  contribution  of  fluctuating  velocity  components  normalized  by  Ufj ,  in  Sy  slices 
at  y/W  =  0  for  AOI  =  30°  and  AOI  =  45°. 
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(f)  AOI  =  15°,  \  v,2/Ujj 
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(g)  AOI  =  0°,  |  w'Z/Ujj 
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(h)  AOI  =  15°,  |  wl2/Ujj 


Figure  4.11:  TKE/f  and  the  contribution  of  fluctuating  velocity  components  normalized  by  U}[ ,  in  Sz  slices 
at  z/H  =  0.25  for  AOI  =  0°  and  AOI  =15°. 
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(e)  AOI  =  30°,  |  v'2/Ujj 


-0.08 

0.06 

0.04 

0.02 


0  12  3  4 

x/S 

(f)  AOI  =  45°,  | 


0.08 

0.06 

0.04 

0.02 

0 


1.5r 


(g)  AOI  =  30°,  \  w'2/U2H 
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Figure  4.12:  TKE#  and  the  contribution  of  fluctuating  velocity  components  normalized  by  Ujj,  in  Sz  slices 
at  z/H  =  0.25  for  AOI  =  0°  and  AOI  =  15°. 
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Figure  4.13:  —u'w ' normalized  by  Ujj,  in  the  Sy  slice  at  y/W  =  0  for  AOI  =  0°,  15°,  30°  and  45°. 
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4.2.4  Arch  Vortex 


Sousa  [2002]  used  the  normalized  angular  momentum  to  locate  the  core  of  the  vortical  structures  behind  a 
single  cuboid  obstacle.  This  normalized  angular  momentum  is  a  3D  implementation  of  Ti  as  defined  by 
Graftieaux  et  al.  [2001].  Monnier  et  al.  [2010]  applied  this  method  to  the  3D  data  of  SPIV  in  street  canyons 
of  an  urban  array  model.  Unlike  the  majority  of  coherent  structures  detection  methods  which  use  the  velocity 
gradients  in  calculations,  Ti  is  an  integral  based  quantity.  T i  is  calculated  as: 


rx 


{XF)  =  hLn 


(; X  —  Xp )  X  V  ( X )  A 

|*-*P||V(*)| 


(4.2) 


where  xp  is  the  spatial  location  at  which  Ti  is  computed,  U  is  the  spatial  domain  over  which  the  integration 
is  performed  (typically  a  small  subset  of  the  entire  data  domain  centered  about  xp)  and  V(x)  is  the  velocity 
vector  at  x  .  The  quantity  Ti  is  therefore  a  vector  quantity. 

We  are  using  the  norm  of  Ti  to  identify  coherent  structures  in  the  urban  array.  The  thresholds  used  to 
display  isosurfaces  of  Ti  are  selected  as  Ti  =  0.4  for  AOI  =  0°  and  15°,  Ti  =  0.35  for  AOI  =  30°,  and 
T i  =  0.33  for  AOI  =  45°.  As  shown  in  Figure  4.15,  the  arch  vortex  (also  known  as  portal  vortex,  Kim 
and  Baik  [2004])  is  easily  observed  in  all  streets  for  AOI  =  0°.  A  reasonable  symmetry  with  respect  to 
the  mid-span  xz  plane  at  y/W  =  0  is  observed,  which  is  expected  for  a  symmetrical  geometry.  This  was 
investigated  in  several  studies  such  as  Monnier  et  al.  [2010],  Becker  et  al.  [2002]  and  Kim  and  Baik  [2004], 
employing  different  geometries.  The  latter  study  used  the  most  similar  geometry  to  the  current  work,  with 
the  geometrical  ratios  of  H/W  =  1  and  S/W  =  1  (compared  with  the  present  ratios  of  H/W  =  2  and 
S/W  =  1.5)  but  the  results  show  a  point  of  contrast  in  general  inclination  of  the  structures. 

For  AOI  =15°  (Figure  4.16)  the  western  leg  in  all  four  streets  moves  toward  the  northwestern  corner 
of  the  street  (moves  closer  to  the  blocks  surface  as  compared  with  AOI  =  0°)  and  the  eastern  leg  moves 
downstream  in  the  street  without  considerable  transition  in  span.  However,  for  the  AOI  =  30°,  results 
presented  in  Figure  4.17,  the  eastern  leg  in  the  1st  and  2nd  streets  is  located  close  to  mid-span  in  the  streets 
and  the  western  leg  is  not  entirely  detectable  within  the  street  and  intersection  boundaries  for  the  thresholds 
used.  Also  the  leg  may  be  outside  of  our  data  acquisition  region.  In  the  3rd  and  4th  streets  the  entire 
structure  is  detected  within  the  street  region;  however  for  the  threshold  used  a  continuous  arch  vortex  in  the 
upper  layer  of  the  3rd  and  4th  streets  is  not  observed.  For  AOI  =  45°  (Figure  4.18)  the  whole  structure  has 
a  transition  upstream  in  the  streets  as  compared  to  AOI  =  30°  results.  In  all  four  streets  only  one  leg  of  the 
arch  (the  eastern  leg)  is  detectable  within  the  street  and  intersection  areas  investigated  in  the  study. 
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(a)  Isometric  view 


(b)  Side  view 


(c)  Top  view 


Figure  4.15:  Arch  vortex,  using  Ti  Iso-surfaces  for  AOI  =  0°. 
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(a)  Isometric  view 


(b)  Side  view 

R£>  ■£■2$  ■ 


(c)  Top  view 


Figure  4.16:  Arch  vortex,  using  Ti  Iso-surfaces  for  AOI  =  15°. 
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(b)  Side  view 


(c)  Top  view 


Figure  4.17:  Arch  vortex,  using  Ti  Iso-surfaces  for  AOI  =  30°. 
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(a)  Isometric  view 
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(b)  Side  view 


(c)  Top  view 


Figure  4.18:  Arch  vortex,  using  Ti  Iso-surfaces  for  AOI  =  45°. 


Chapter  5 


Numerical  results 


We  start  this  chapter  with  a  discussion  of  several  code  validation  and  test  cases  demonstrating  accuracy  and 
performance  of  our  computer  code.  After  that,  the  following  sections  describe  simulation  results  for  the 
urban  boundary  layer  flow. 

5.1  Code  validation 

5.1.1  Wannier  Flow 

Our  first  test  case  is  the  Wannier  flow,  a  Stokes  flow  past  a  rotating  cylinder  next  to  a  sliding  wall.  The 
schematic  of  the  flow  domain  and  streamlines  are  plotted  in  Figure  5.1.  The  solution  can  be  described  in 
terms  of  the  cylinder  radius  (r),  its  angular  velocity  (cc),  distance  from  the  center  of  the  cylinder  to  the 
moving  wall  ( d )  and  the  velocity  of  the  wall  ( U ).  In  this  simulation  a  cylinder  of  radius  r  =  0.25  is  centered 
at  (xi,^2)  =  (0,0)  and  is  rotating  with  an  angular  velocity  cc  =  2.  The  wall  is  located  at  d  =  0.5  and 
moves  with  velocity  U  =  1.0.  The  exact  solution,  originally  derived  by  Wannier  [1950]  and  Karniadakis 
and  Sherwin  [2005],  can  be  written  as: 


(5.1) 


(5.2) 
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where 


s  =d2  —  r2, 

T 

u 

a°  =  ln(r)  ’ 

f  1  r2w\ 

a  i 

<22  =2(<i  +  s)  | 

vao+2~J’ 

03 
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Figure  5.1:  Schematic  of  Wannier  flow. 


Figure  5.2:  Mesh  with  201  quadrilateral  spectral  elements. 

The  mesh  used  for  the  simulation  is  shown  in  Figure  5.2.  The  outer  mesh  is  made  of  201  quadrilateral 
spectral  elements  and  the  polynomial  order  is  varied  from  4  to  10  for  this  study. 
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Figure  5.3  plots  the  and  L 2  norms  of  error  in  velocity.  The  L2  norms  of  error  are  normalized  by  the 
I/2  norms  of  the  exact  solution  for  each  of  the  velocity  components  to  obtain  relative  errors.  We  can  see  that 
the  errors  in  both  norms  drop  exponentially  fast.  At  higher  values  of  polynomial  degree,  the  errors  are 
dominated  by  errors  at  points  close  to  the  cylinder  wall  whereas  the  L2  error  continues  to  drop  exponentially. 
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Figure  5.3:  Convergence  in  and  L2  norm  as  a  function  polynomial  order  for 
Wannier  flow. 
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5.1.2  Kovasznay  Flow 

Our  second  test  case  is  Kovasznay  flow,  a  laminar,  incompressible  flow  behind  a  two-dimensional  grid.  The 
exact  solution,  due  to  Kovasznay  [1948],  can  be  written  as  a  function  of  Reynolds  number  (Re)  as  follows: 


where 


u(x,y)  = 

1  —  eXxcos(2ny), 

(5.3) 

v(x,y)  = 

2^eXxsin(2ny), 

Z7T 

(5.4) 

Re 

(Re2  2V 

a  =  T“ 

br+47r  • 

(5.5) 

The  schematic  of  the  flow  domain  is  shown  in  Figure  5.4.  Dirichlet  boundary  conditions,  obtained  from 
the  exact  solution,  are  prescribed  on  all  domain  boundaries.  Second-order  time  stepping  is  used  and  the 
solution  is  marched  until  it  reaches  a  steady  state.  Various  parameters  used  for  the  simulation  are  presented 
in  Table  5.1.  Figure  5.5  shows  a  typical  mesh  used  for  the  simulation  of  Kovasznay  flow.  Figure  5.6  shows 
the  steady-state  streamlines  obtained  from  the  simulation.  In  Figure  5.7  we  see  that  the  L0 c  and  L 2  norm  of 
error  both  decrease  exponentially  as  we  increase  the  order  of  polynomial  expansion  (P).  As  in  the  earlier 
case,  I/2  norms  of  error  are  normalized  with  the  L2  norm  of  the  exact  solution. 


Periodic 


Periodic 


Figure  5.4:  Schematic  of  Kovasznay  flow. 

Figure  5.8  plots  the  effect  of  stream  wise  domain  size  on  solution  accuracy  when  outflow  boundary  condi¬ 
tions  are  used  at  the  right  boundary.  Initially,  we  can  see  that  spectral  convergence  is  obtained  with  increase 
in  polynomial  degree.  After  a  certain  point  the  error  saturates  as  the  error  due  to  outflow  boundary  condi¬ 
tion  dominates  numerical  errors.  This  can  be  delayed  by  moving  the  outflow  boundary  further  downstream. 
The  simulation  is  also  performed  using  a  3D  mesh  shown  in  Figure  5.9.  Periodic  boundary  conditions  are 
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Table  5.1:  Parameters  for  Kovasznay  flow  simulation. 


Parameter 

Value 

Re 

40.0 

At 

0.001 

Ncells 

12 

P 

4-14 

border 

2 

used  in  the  third  direction.  As  in  the  2D  case,  Figure  5.10  shows  exponential  convergence  with  increase  in 
polynomial  order. 
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Figure  5.5:  Mesh  used  for  the  simulation  of  steady  state  Kovasznay  flow 
(Ncells  =  12;  P=14.). 
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Figure  5.6:  Streamlines  for  steady  state  Kovasznay  flow. 
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Figure  5.7:  Convergence  in  and  L 2  norm  as  a  function  polynomial  order  for  steady 
state  Kovasznay  flow  at  Reynolds  number  40. 


(a)  (b) 

Figure  5.8:  Study  of  effect  of  domain  size  on  solution  accuracy  for  2d  Kovasznay  flow 
with  outflow  boundary  condition  prescribed  on  the  right  face,  (a)  L ^  and  L2 
errors  in  u  and  v  as  a  function  of  polynomial  order  for  domain  x  G  [—0.5  5.0]. 

(b)  I/2  errors  in  u  as  a  function  of  polynomial  order  for  various  domain  sizes. 
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Figure  5.9:  3-d  mesh  used  for  Kovasznay  flow  simulation  at  Reynolds  number  40. 
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Figure  5.10:  Convergence  in  and  L2  norm  as  a  function  polynomial  order  for  3-d 
steady  state  Kovasznay  flow  at  Reynolds  number  40. 
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5.1.3  Orr-Sommerfeld  problem 

Our  third  test  case  is  the  simulation  of  the  highly  sensitive  Orr-Sommerfeld  problem.  This  involves  compu¬ 
tation  of  growth  rates  of  small-amplitude  Tollmien-Schlichting  waves  in  the  plane  Poiseuille  flow  (Fischer 
[1997]).  For  this  case,  an  accurate  solution  is  available  from  linear  stability  theory.  The  domain  and  the 
mesh  used  for  this  simulation  is  shown  in  Figure  5.11.  The  domain  consists  of  two  walls  separated  by  a  dis¬ 
tance  2h  =  2  and  a  stream  wise  length  of  27 rh.  Periodic  boundary  conditions  are  used  in  the  stream  wise  di¬ 
rection  and  the  flow  is  driven  by  a  constant  body  force.  The  Reynolds  number  is  given  by  Re  =  =  7500. 

Also  shown  is  the  perturbation  stream  function.  The  initial  condition  is  the  solution  for  the  plane  Poiseuille 
flow  superimposed  with  a  perturbation.  Specifically,  the  initial  condition  is  given  by 

u  =  1  —  y2  +  eft,  (5.6) 

v  =  ev.  (5.7) 


Here  (u,  v)  correspond  to  the  only  unstable  eigenfunction  of  the  Orr-Sommerfeld  equation  with  wave 
number  unity  at  Re  =  7500.  We  use  e  =  10-5. 

Linear  stability  theory  predicts  the  energy  of  the  perturbation 

p2ir  p  1 

E(t)  =  J  J  [(1 -y2 -u)2 +v2]dydx,  (5.8) 


to  grow  as  e2wt,  where  c u  =  0.002234975649.  The  relative  error  in  growth  rate  is  given  by 


error  = 


1 

UJ 


{  E(t  +  to)  \ 

V  E(to)  J 


(5.9) 


Table  5.2  shows  the  computed  perturbation  energy  for  various  polynomial  degrees.  Figure  5.12  plots  the 
temporal  evolution  of  the  perturbation  and  the  computed  growth  rate  for  P  =  9.  In  all  cases,  to  is  chosen 
to  be  1  sec  after  the  start  of  the  simulation.  A  second-order  accurate  time  stepping  scheme  with  a  time  step 
of  At  =  0.003125  is  used  for  these  simulations.  For  each  of  the  cases,  the  relative  error  in  growth  rate  is 
computed  after  t  =  50.28  sec  which  corresponds  to  two  periods  of  oscillation  for  the  Tollmien-Schlichting 
waves.  We  can  see  that  the  computed  value  of  the  growth  rate  is  in  excellent  agreement  with  the  result  from 
linear  stability  theory.  We  can  also  see  spectral  convergence  as  the  polynomial  degree  is  increased  from  7 
to  13.  Once  the  polynomial  degree  is  higher  than  13,  the  temporal  error  corresponding  to  the  time  step  size 
dominates  the  spatial  error. 


Figure  5.11:  Domain  and  spectral  element  mesh  used  for  the  Orr-Sommerfeld 
simulation.  Also  shown  is  the  stream  function  of  the  perturbation. 


69 


2.238 


x  10"' 


q:poooooaxxxxxrx^^ 


2.236 


2.234 


3  2.232 


2.23 


2.228 

C) 

2.226  — 


100 


200  300 

time  (sec) 


400 


500 


Figure  5.12:  Temporal  evolution  of  perturbation  energy  (left)  and  growth  rate  (right) 
for  P=9,  At  =  0.003125. 


Table  5.2:  Spatial  convergence,  Orr-Sommerfeld  problem:  Nei  =  15, 
At  =  0.003125. 


Degree 

UJ 

error 

7 

2.293545 1 88889385e-03 

2.6205 8962099 1 660e-02 

9 

2.237078145708761e-03 

9.407246605581612e-04 

11 

2.234894879598453e-03 

3.613882844038103e-05 

13 

2.234950158073417e-03 

1 . 1 4054605449 1 922e-05 

15 

2.234950158073417e-03 

1 . 1 4054605449 1 922e-05 
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5.1.4  Vortex  shedding  from  circular  cylinder 

Our  next  test  case  is  laminar  flow  past  a  circular  cylinder  in  free  stream.  For  Re  >  47,  vortex  shedding  is 
known  to  occur  and  a  Von-Karman  vortex  street  forms  in  the  wake  of  the  cylinder.  This  vortex  shedding 
causes  oscillating  forces  on  the  cylinder  at  a  particular  frequency,  called  the  Strouhal  frequency,  and  is  of 
particular  importance  in  engineering  design.  In  our  current  study,  we  perform  this  experiment  at  various 
Reynolds  numbers  and  compare  the  results  with  other  simulations  and  experimental  data.  Strouhal  number 
is  used  as  a  measure  to  compare  the  data. 

The  domain  used  for  this  simulation  is  shown  in  Figure  5.13.  Also  shown  is  the  spectral  element 
mesh  used  for  the  simulation  at  Re  =  100.  The  mesh  is  further  refined  for  higher  Reynolds  number 
cases.  The  diameter  of  the  cylinder  (d)  is  chosen  as  d  m  0.2828.  The  center  of  the  cylinder  is  located 
at  (#1,  £2)  =  (1-5,  0.5).  The  domain  starts  at  about  5  diameters  upstream  from  the  center  of  the  cylinder 
and  extends  up  to  16  diameters  downstream  from  the  center  of  the  cylinder.  Uniform  inflow,  =  1.0,  is 
prescribed  at  the  inflow  boundary  and  an  outflow  boundary  condition  is  prescribed  at  the  outflow  boundary. 
Side  boundaries  are  located  at  about  5  diameters  from  the  center  of  the  cylinder.  A  symmetry  boundary 
condition  is  prescribed  on  these  boundaries.  The  Reynolds  number  ,  Re  =  U^d/v ,  is  varied  from  47  to 
387.  Figure  5.14  plots  the  instantaneous  velocity  and  pressure  contours  at  Re  =  100. 

Figure  5.15  compares  the  results  of  the  simulation  with  2D  direct  numerical  simulations  by  Henderson 
[1997]  and  with  the  experiments  of  Williamson  [1989].  The  Strouhal  numbers  obtained  at  various  Reynolds 
numbers  are  in  good  agreement  with  2D  simulations  of  Henderson  at  all  Reynolds  numbers  considered  in 
this  study.  They  are  also  in  good  agreement  with  the  experimental  data  for  Re  <  190  at  which  the  2D 
wake  becomes  unstable  and  bifurcates  to  a  three-dimensional  flow.  Above  this  Reynolds  number,  three 
dimensional  simulations  are  required  to  accurately  resolve  all  flow  features. 


Figure  5.13:  Spectral  element  mesh  for  simulation  of  vortex  shedding  from  a  circular 
cylinder  at  Re=100. 
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(c)  Pressure 


Figure  5.14:  Instantaneous  velocity  and  pressure  contours  at  Re  =  100. 


Figure  5.15:  Comparison  of  variation  of  Strouhal  number  with  Reynolds  number. 
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5.1.5  Flow  over  a  backward  facing  step 

Our  final  test  case  is  the  three-dimensional  laminar  flow  over  a  backward-facing  step  with  an  expansion  ratio 
of  1:1.94.  This  is  based  on  the  experimental  setup  of  Armaly  et  al.  [1983]  and  was  previously  used  in  the 
numerical  simulations  of  Couzy  [1995],  Biswas  et  al.  [2004]  and  Shahbazi  [2007].  The  flow  is  characterized 
by  three  recirculation  regions  whose  occurrence,  location  and  size  depends  on  the  Reynolds  number  and 
expansion  ratio  of  the  step.  In  this  study,  we  restrict  ourselves  to  Re  <  400  which  has  a  single  recirculation 
zone  on  the  bottom  wall  directly  at  the  base  of  the  step. 

The  geometry  and  mesh  used  for  the  simulation  is  shown  in  Figure  5.16.  The  height  of  the  downstream 
section,  h ,  is  chosen  to  be  1.  The  inflow  is  located  3 h  upstream  of  the  expansion  and  the  outflow  is  lo¬ 
cated  20 h  downstream  of  the  expansion.  The  spanwise  length  is  9 h  with  a  wall  on  one  side  and  symmetry 
boundary  conditions  on  the  opposite  side.  The  inlet  velocity  profile  is  chosen  to  be  the  tensor  product  of  a 
parabola  (in  the  z-direction)  and  a  Blasius  boundary  layer  (in  the  ^/-direction).  The  Blasius  velocity  profile 
is  characterized  by  the  boundary  layer  thickness  (5)  which  represents  the  wall-normal  distance  at  which  the 
velocity  attains  99%  of  the  free-stream  value.  Specifically,  the  inlet  velocity  is  given  by 


u  =  [15.08739(0.5149  -  z)z\  b(y ), 

(5.10) 

v  =  0.0, 

(5.11) 

w  =  0.0. 

(5.12) 

The  mesh  consists  of  682  spectral  elements  with  with  ninth-degree  polynomial  approximation  within  each 
element.  Second-order  time  stepping  with  A  =  0.0025  is  used.  The  simulation  is  continued  until  a  steady- 
state  solution  is  obtained.  The  simulation  is  carried  out  at  Reynolds  numbers  172  and  343.  The  recirculation 
zones  for  these  Reynolds  numbers  are  plotted  in  Figures  5.17  and  5.18.  Figure  5.19  shows  the  variation  of 
the  length  of  the  recirculation  region  with  Reynolds  number.  We  can  see  that  the  length  of  the  recirculation 
zone  increases  with  increase  in  Reynolds  number.  These  results  are  in  excellent  agreement  with  earlier 
numerical  simulations  and  with  the  experimental  data. 


Figure  5.16:  Spectral  element  mesh  for  simulation  of  flow  over  a  backward 
facing  step. 
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Figure  5.17:  Flow  pattern  near  the  symmetry  plane  of  the  backward  facing  step  for 
Re=172. 
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Figure  5.18:  Flow  pattern  near  the  symmetry  plane  of  the  backward  facing  step  for 
Re=343. 
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Figure  5.19:  Variation  of  length  of  primary  recirculation  region  with  Reynolds  number 
for  flow  over  a  backward  facing  step. 
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5.2  The  urban  boundary  layer  experiment 

5.2.1  The  urban  boundary  layer 

The  atmospheric  boundary  layer  (ABL)  is  defined  in  (Stull  [1988])  as  the  part  of  the  troposphere  that  is 
directly  influenced  by  the  presence  of  the  Earth’s  surface  and  responds  to  surface  forcing  with  a  time  scale  of 
about  an  hour  or  less.  It  consists  of  a  succession  of  quasi-equilibrium  boundary  layers  diffusing  within  older 
boundary  layers  in  response  to  the  forcing  by  frictional  drag,  terrain,  heat  transfer  and  pressure  gradients. 
The  lowest  part  of  the  ABL  is  called  the  surface  layer.  It  extends  to  a  height  of  about  150  m  depending  on 
the  terrain.  The  surface  layer  over  an  urban  area  can  be  classified  into  three  sublayers:  the  urban  canopy 
layer,  the  roughness  sublayer  and  the  inertial  sublayer  (Britter  and  Hanna  [2003]).  This  is  schematically 
shown  in  Figure  5.20.  In  the  urban  canopy  layer  the  flow  at  a  specific  point  is  directly  influenced  by  the  local 
obstacles,  and  in  the  roughness  sublayer  the  flow  is  still  adjusting  to  the  effect  of  many  obstacles. 


suTts-c^ 

muchr.es  a 

Eubla/or 


k 


Figure  5.20:  Schematic  of  Flow  through  and  over  an  urban  area  (Grimmond  and  Oke  [1999]). 

The  inertial  sublayer  is  the  area  where  the  boundary  layer  has  adapted  to  the  integrated  effect  of  the 
underlying  surface.  The  mean  velocity  profile  in  this  part  can  be  represented  by  the  log-law: 

“=0),n 


z  —  d 
zo 


(5.13) 
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where  u*  is  the  friction  velocity,  z  is  the  elevation,  k  is  the  von=Karman  constant,  zo  is  the  surface  roughness 
length  and  d  is  the  surface  displacement  length.  The  last  two  parameters  depend  on  the  terrain,  increasing 
with  the  roughness  of  the  terrain.  Various  parameterizations  for  zq  based  on  land-use  were  presented  in  the 
works  by  Stull  [1988],  Theurer  [1999],  Grimmond  and  Oke  [1999]  and  Davenport  [1965].  An  alternative 
representation  of  velocity  profile  in  the  inertial  layer  widely  used  in  wind  profile  modeling  is  through  a 
power  law 


(5.14) 


where  u§  is  the  mean  velocity  at  z  =  5  and  the  exponent  a  depends  on  the  roughness  of  the  surface  (Wang 
et  al.  [1996]).  This  can  be  reformulated  as 


u 


Ui 


(5.15) 


where  Ui  is  the  velocity  at  a  point  y  —  yi  in  the  inertial  layer.  The  advantage  of  this  power  law  is  that 
it  involves  only  one  scaling  parameter.  While  velocity  profiles  in  the  inertial  layer  can  be  satisfactorily 
described  by  either  a  power  law  or  a  log  law  with  appropriately  chosen  parameters,  flows  inside  urban 
canopies  are  in  general  three-dimensional  and  quite  complex  and  are  the  focus  of  our  current  work. 

Figure  5.21  shows  a  schematic  of  flow  around  a  surface-mounted  cube  from  Martinuzzi  and  Tropea 
[1993].  This  is  a  fundamental  building  block  for  understanding  flow  characteristics  in  urban  street  canyons. 
In  this  case,  a  cube  is  placed  in  a  fully-developed  channel  flow.  The  cube  imposes  a  strong  adverse  pressure 
gradient  on  the  flow  causing  boundary  layer  separation  upstream  of  the  cube.  This  causes  the  development 
of  strong  recirculation  regions  upstream  and  around  the  cube.  The  approach  boundary  layer  vorticity  rolls 
up  into  a  horse  shoe  vortex  which,  along  with  the  arch  vortex  located  immediately  behind  the  cube,  is  one  of 
the  most  prominent  features  of  this  flow.  A  favorable  pressure  gradient  diverts  the  flow  upstream  of  the  cube 
from  the  symmetry  plane  towards  the  sides.  Shear  layers  separate  on  top  and  sides  of  the  block.  Behind 
the  cube,  there  is  a  favorable  pressure  gradient  from  the  sides  to  the  symmetry  plane  causing  the  flow  to 
reattach  behind  the  block.  Sousa  [2002]  studied  the  turbulent  flow  around  a  surface-mounted  obstacle  using 
two-dimensional  three-component  DPIV  (Digital  Particle  Image  Velocimetry).  The  out-of-plane  velocity 
component  was  obtained  by  the  use  of  continuity  applied  to  two-dimensional  velocity  fields  recorded  in 
parallel  planes.  He  noted  that  the  use  of  swirling  strength  and  normalized  angular  momentum  to  identify 
vortices  is  superior  to  traditional  vorticity-based  methods. 

Becker  et  al.  [2002]  investigated  the  effect  of  wind  direction,  aspect  ratio,  Reynolds  number,  and  the 
boundary  layer  type  on  flow  structures  around  a  single  obstacle.  This  work  extends  the  arch  vortex  topology 
proposed  by  Martinuzzi  and  Tropea  [1993]  to  non-zero  angle  of  attack.  They  noticed  that  increasing  the 
angle  of  attack  caused  a  dislocation  of  one  of  the  footprints  of  the  vortex  until  it  switched  to  the  top  of  the 
obstacle  at  an  angle  of  60° .  They  found  no  fundamental  difference  in  the  vortex  structure  for  various  bound¬ 
ary  layers  and  found  that  the  length  of  the  recirculation  region  decreased  for  rougher  incoming  boundary 
layers. 

Additional  complexities  arise  when  an  object  is  placed  in  the  wake  of  another  object.  Based  on  the 
reattachment  length  of  the  wake  behind  the  isolated  object  (in  the  absence  of  objects  downstream)  and  the 
downstream  spacing  of  the  second  object,  the  flow  can  be  broadly  classified  into  three  different  regimes. 
Figure  5.22  shows  the  classification  into  isolated  roughness,  wake  interference  and  skimming  flow  regimes 
based  on  the  canyon  aspect  ratio  (Oke  [1988]).  Here,  we  use  the  words  “canyon”  and  “streets”  interchange¬ 
ably.  The  canyon  aspect  ratio  is  defined  as  the  ratio  of  spacing  between  the  blocks  to  the  height  of  the  block 
(j|).  This  work  is  2D  in  nature  and  assumes  infinite  width  of  the  street  canyon.  When  the  canyon  aspect 
ratio  is  small  (j|  <  1.5),  the  majority  of  the  flow  skims  over  the  canyon  with  vortices  trapped  within  the 
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canyon.  This  is  called  the  skimming  flow  regime.  When  the  objects  are  spaced  further  apart  but  not  too  far 
(1.5  <  <  8  —  9),  the  downstream  object  “sees”  the  wake  of  the  upstream  object.  This  is  called  the  wake 

interference  regime.  Finally,  when  the  canyon  aspect  ratio  is  larger  d  >  8  —  9)  the  individual  elements 
act  as  isolated  roughness  elements  and  the  interaction  of  the  flows  induced  by  the  buildings  is  negligible. 
This  is  called  the  isolated  roughness  regime.  ?  have  studied  the  effect  of  street  canyon  aspect  ratio  on  flow 
characteristics  using  a  water  tunnel  experiment.  Laser  Doppler  anemometry  is  used  for  measuring  velocity. 
They  noticed  that  similar  regimes  were  found  even  when  the  width  of  the  street  is  finite. 


Figure  5.21:  Schematic  of  Flow  around  a  surface  mounted  cube  (Martinuzzi  &  Tropea  [1993]). 


(a)  Isolated  roughness  flow 


(b)  Wake  interference  flow 


(c )  Skimming  flow 


Figure  5.22:  Schematic  illustrating  (a)  isolated  roughness  (b)  wake  interference  and  (c)  skimming  flow 
regimes  in  an  urban  street  canyon  by  Oke  [1988]. 


The  MUST  (Mock  Urban  Setting  Test)  experiment  (Biltoft  [2001])  is  a  scaled  urban  dispersion  experi¬ 
ment  performed  at  Dugway  proving  grounds,  Utah  in  2001.  The  objective  was  to  overcome  the  scaling  and 
measurement  limitations  of  laboratory  experiments  and  characterize  difficulties  presented  by  real  urban  set- 
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tings.  A  10  x  12  array  of  containers  was  used  to  represent  an  urban  domain  with  symmetric  characteristics. 
Continuous  tracer  releases  were  generated  with  propylene  as  a  tracer  gas,  and  concentrations  were  measured 
using  fast  response  photo-ionization  detectors.  Gailis  [2004]  conducted  a  wind  tunnel  dispersion  study  of 
the  MUST  array.  This  study  was  intended  to  bridge  some  of  the  gaps  between  laboratory  and  full-scale  out¬ 
door  trials.  Yee  et  al.  [2006]  compared  the  field  experiment  to  both  a  l:50-scale  wind  tunnel  and  l:205-scale 
water  channel  simulation.  The  study  notes  that  while  the  behavior  of  plume  dispersion  is  qualitatively  the 
same  in  all  studies,  appropriately  scaled  water-channel  simulations  were  able  to  reproduce  qualitatively  the 
results  of  the  full-scale  field  experiments  better  than  wind-tunnel  experiments. 

Monnier  et  al.  [2010]  has  investigated  the  three-dimensional  flow  through  an  urban-type  array  consisting 
of  rows  of  cuboid  plexiglas  blocks  in  a  laboratory-modeled  neutrally  stratified  atmospheric  boundary  layer 
using  SPIV  (stereoscopic  particle  image  velocimetry).  A  typical  domain  used  in  the  experiment  is  shown  in 
Figure  5.23.  The  domain  can  be  defined  in  terms  of  width  of  the  block  (TV),  height  of  the  block  ( H )  and 
spacing  between  the  blocks  ( S ).  The  SPIV  technique  is  used  to  acquire  velocity  data  in  vertical  planes  in 
between  the  blocks.  In  this  experiment,  the  effect  of  stream- wise  spacing  between  adjacent  rows  defining 
the  two  different  flow  regimes  (wake  interference  and  skimming  flow  regimes)  as  well  as  the  effect  of 
the  incident  angle  of  the  approaching  boundary  layer  are  studied.  Dominant  mechanisms  responsible  for 
transport  and  dispersion  were  quantified. 

While  earlier  work  in  this  area  is  primarily  based  on  field  measurements  and  wind  tunnel  experiments, 
rapid  increase  in  computing  power  and  development  of  highly  scalable  parallel  algorithms  to  harness  this 
power  is  making  numerical  simulations  a  tool  of  choice.  Earliest  numerical  simulations  in  this  area  modeled 
flow  in  an  archetypal  street  canyon  which  is  basically  a  turbulent  shear  flow  above  a  rectangular  cavity  with 
mean  flow  perpendicular  to  the  axis  of  the  street  canyon.  While  2D  simulations  assumed  infinite  street  width, 
3D  simulations  used  periodic  boundary  conditions  in  the  span  wise  direction.  Baik  and  Kim  [1999]  modeled 
flow  and  pollution  in  a  2D  street  canyon  using  a  2D  ft  —  e  turbulence  model .  Their  code  used  a  finite  volume 
method  with  a  staggered  grid.  A  power-law  velocity  profile  is  used  as  an  inlet  boundary  condition  at  the  top 
of  the  upstream  building.  They  studied  the  flow  patterns  in  a  street  canyon  for  various  street  canyon  aspect 
ratios.  Kim  and  Baik  [2004]  performed  3D  numerical  simulations  of  flow  within  an  array  of  cube  using  a 
renormalization  group  k  —  e  scheme.  They  investigated  the  effect  of  angle  of  incidence  on  flow  structures 
and  classified  the  flow  in  three  regimes  based  on  the  angle  of  incidence.  Liu  et  al.  [2004]  performed  large 
eddy  simulations  of  flow  in  a  model  urban  street  canyon  using  the  Smagorinsky  subgrid  scale  model.  They 
investigated  the  effect  of  street  canyon  aspect  ratio  on  flow  structures  and  pollutant  transport  within  the  street 
canyon.  Shah  and  Ferziger  [1997]  studied  flow  over  a  surface  mounted  cube  using  a  large  eddy  simulation.  A 
second-order  finite  volume  code  was  used.  Data  from  large-eddy  simulations  of  channel  flow  at  comparable 
Reynolds  number  is  used  to  generate  boundary  conditions  upstream  of  the  cube.  Camelli  and  Lohner  [2006] 
studied  flow  and  dispersion  patterns  in  realistic  urban  areas  like  the  Tyson’s  corner  area  in  Fairfax,  and  the 
Madison  Square  Garden  area  in  New  York  City  using  VLES  (Very  Large  Eddy  Simulation).  They  used  a 
finite  element  code  with  a  dynamic  Smagorinsky  LES  model.  Tseng  et  al.  [2006]  studied  flow  and  dispersion 
through  a  model  of  downtown  using  large  eddy  simulation.  Their  code  uses  a  pseudo- spectral  method  in  the 
horizontal  directions  and  a  second-order  accurate  central  difference  scheme  in  the  vertical  direction.  The 
presence  of  bluff  bodies  is  modeled  using  an  immersed  boundary  method.  A  Lagrangian  dynamic  LES 
(large  eddy  simulation)  model  (Bou-Zeid  et  al.  [2005])  is  used  for  modeling  sub-grid  scale  stresses. 

In  this  section,  we  present  the  results  of  the  spectral  element  simulation  of  flow  in  a  model  urban  street 
canyon.  The  domain  used  in  this  simulation  consists  of  a  5  x  7  array  of  blocks  and  is  similar  to  that  used  in 
the  work  of  Monnier  et  al.  [2010].  The  inflow  velocity  profile  is  prescribed  based  on  data  provided  from  the 
hot-wire  measurements  and  the  flow  upstream  of  the  array  is  accurately  resolved.  Low  pass  filtering  is  used 
for  stabilizing  the  simulation.  No  LES  model  is  used. 
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In  section  5.2.2  the  results  for  the  0°  AOI  (angle  of  incidence)  case,  where  the  approaching  boundary 
layer  is  normal  to  the  array  of  blocks,  are  presented  in  detail.  These  results  are  compared  with  experimental 
data.  Some  preliminary  results  for  15°  AOI  are  presented  in  section  5.2.3. 


Figure  5.23:  Array  of  cuboid  plexiglas  blocks  used  in  the  experiment  by 
Monnier  et  al.  [2010]. 


5.2.2  Numerical  simulation  of  the  urban  boundary  layer 
at  zero  angle  of  incidence 

Description  of  numerical  simulation  and  comparison  with  experimental  data 

The  model  of  the  urban  street  canyon  used  in  the  current  study  is  shown  in  figure  5.24.  It  consists  of  a 
5x7  array  of  blocks.  Each  block  is  1  unit  long  ( L  =  1),  1  unit  wide  (W  =  1)  and  2  units  high  (H  =  2), 
respectively.  The  spanwise  and  streamwise  spacing  between  blocks  is  1.5  units  ( S  =  1.5).  The  wind 
direction  is  normal  to  the  array  ( AOI  =  0°).  The  domain  used  for  our  numerical  simulation  is  shown  in 
figure  5.25.  The  domain  consists  of  a  single  row  of  blocks  and  periodic  boundary  conditions  prescribed 
in  the  spanwise  direction.The  Reynolds  number  ( Ren ),  based  on  block  height  and  inlet  velocity  at  block 
height,  is  Ren  =  UhH/u  =  6283.  The  domain  used  for  the  simulation  is  shown  in  figure  5.25.  The  domain 
consists  of  a  single  row  of  blocks  and  periodic  boundary  conditions  are  used  in  the  spanwise  direction.  The 
inflow  boundary  is  located  10  block  widths  upstream  of  the  first  block  and  the  outflow  boundary  is  located 
40  block  widths  downstream  of  the  last  block.  A  symmetry  boundary  is  condition  is  used  on  the  top  and 
is  located  at  17  block  widths  above  the  bottom  wall.  Inlet  velocity  profiles  obtained  from  the  hot-wire  data 
provide  the  inlet  boundary  conditions  for  the  simulation. 

Three  different  simulations  were  carried  out  with  the  objective  of  studying  the  effects  of  grid  resolution. 

The  first  simulation  uses  a  mesh  of  15488  hexahedral  spectral  elements  with  a  degree-6  polynomial 
approximation  within  each  spectral  element.  This  corresponds  to  about  5.3  million  degrees  of  freedom.  This 
mesh  is  shown  in  Figure  5.26.  Third-order  time  stepping  is  used.  The  backward  difference  scheme  is  used 
for  diffusion  terms  and  extrapolation  is  used  for  nonlinear  terms.  De-aliasing  and  low  pass  filtering  is  used 
for  stabilizing  the  simulation.  No  turbulence  modeling  is  used.  We  have  used  the  Nek5000  code  for  this 
simulation. 
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Figure  5.24:  A  typical  urban  street  canyon. 


Figure  5.25:  Domain  used  for  the  simulation  of  the  urban  boundary  layer  experiment. 


The  second  simulation  uses  a  mesh  of  45748  hexahedral  spectral  elements  with  a  degree-6  polynomial 
approximation  within  each  spectral  element.  This  corresponds  to  about  15.6  million  degrees  of  freedom. 
This  mesh  is  shown  in  Figure  5.27.  Second-order  time  stepping  is  used.  The  backward  difference  scheme  is 
used  for  diffusion  terms  and  extrapolation  is  used  for  nonlinear  terms.  De-aliasing  and  low  pass  filtering  is 
used  for  stabilizing  the  simulation.  Steady  inflow  boundary  condition  is  used  for  this  case.  This  simulation 
is  performed  using  Spec  solve. 

The  third  simulation  uses  a  mesh  of  88688  spectral  elements  with  a  degree-6  approximation  within  each 
spectral  element.  This  mesh  is  shown  in  Figure  5.28.  This  corresponds  to  about  30.4  million  degrees  of 
freedom. 
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Figure  5.26:  15488  element  hexahedral  mesh  for  simulations  1. 
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(c)  Closeup  of  mesh  along  y=1.25  plane. 


Figure  5.27:  45748  element  hexahedral  mesh  for  simulation  2. 
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(b)  Slice  of  mesh  along  z=0.5.  plane 


(c)  Closeup  of  mesh  along  y=1.25.  plane 


(d)  Closeup  of  mesh  along  z=0.5.  plane 


Figure  5.28:  88688  element  hexahedral  mesh  for  simulation  3. 
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We  compare  the  simulation  results  with  PIV  data  obtained  from  the  experiment  in  street  2.  Comparisons 
are  made  along  the  planes  x  =  0.75,  y  =  1.25  and  2  =  0.5.  Figures  5.29,  5.30  and  5.31  compare  the 
mean  velocity,  rms  velocity  and  Reynolds  stress  between  the  simulation  and  PIV  data  in  the  vertical  plane 
x  =  0.75.  We  can  see  that  there  is  good  agreement  between  simulation  data  and  PIV  data.  Figures  5.32, 
5.33  and  5.34  compare  the  mean  velocity,  rms  velocity  and  Reynolds  stress  between  the  simulation  and 
PIV  data  in  the  vertical  plane  y  =  1.25.  In  this  case,  one  can  notice  a  significant  improvement  in  W  and 
TKE  prediction  between  simulation  1  and  simulation  2.  This  tells  us  that  the  mesh  used  for  simulation  1 
is  not  fine  enough  to  resolve  the  flow  features.  There  is  only  little  improvement  between  simulation  3 
and  simulation  2.  In  general,  all  simulations  seem  to  over-predict  the  magnitude  of  W  upstream  of  the 
third  block.  Figures  5.35,  5.36  and  5.37  compare  the  mean  velocity,  rms  velocity  and  Reynolds  stress 
between  the  simulation  and  PIV  data  on  plane  z  =  0.5.  Figure  5.38  plots  the  first  4  POD  (proper  orthogonal 
decomposition  modes)  for  U,  V  and  W.  Overall,  we  can  see  that  there  is  good  agreement  between  simulation 
and  PIV  data.  From  the  data  obtained  from  the  grid  resolution  study,  we  can  see  that  meshes  two  and  three 
are  fine  enough  to  resolve  all  flow  features.  We  think  the  comparison  can  be  further  improved  by  using  more 
accurate  inflow  boundary  conditions.  Finally,  Figure  5.39  presents  the  scaling  results  for  simulation  3.  We 
notice  the  parallel  efficiency  of  71.46%  as  the  number  of  processors  is  increased  from  252  (  120710  degrees 
of  freedom  per  processor)  to  504(  60360  degrees  of  freedom  per  processor).  These  tests  are  performed  on  the 
Kraken  supercomputer  at  NICS  (National  institute  for  computational  sciences)  through  an  XSEDE  startup 
allocation. 
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Figure  5.29:  Comparison  of  mean  velocity  contours  between  spectral  element 
simulations  and  PIV  data  on  vertical  plane  X  =  0.75. 
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PIV  Simulation  3  Simulation  2  Simulation  1 


Figure  5.30:  Comparison  of  root  mean  square  velocity  and  turbulent  kinetic  energy 
contours  between  spectral  element  simulations  and  PIV  data  on  vertical  plane 


X  =  0.75. 
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Figure  5.31:  Comparison  of  Reynolds  stress  data  between  spectral  element  simulations 
and  PIV  data  on  vertical  plane  X  =  0.75. 
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Figure  5.32:  Comparison  of  mean  velocity  contours  between  spectral  element 
simulations  and  PIV  data  on  vertical  plane  Y  =  1.25. 


91 


PIV  Simulation  3  Simulation  2  Simulation  1 
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Figure  5.33:  Comparison  of  root  mean  square  velocity  and  turbulent  kinetic  energy 
contours  between  spectral  element  simulations  and  PIV  data  on  vertical  plane 
Y  =  1.25. 
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Figure  5.34:  Comparison  of  Reynolds  stress  data  between  spectral  element  simulations 
and  PIV  data  on  vertical  plane  Y  =  1.25. 
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Figure  5.35:  Comparison  of  mean  velocity  contours  between  spectral  element 
simulations  and  PIV  data  on  vertical  plane  Z  =  0.5. 
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Figure  5.36:  Comparison  of  root  mean  square  velocity  and  turbulent  kinetic  energy 
contours  between  spectral  element  simulations  and  PIV  data  on  vertical  plane 


Z  =  0.5. 
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PIV  Simulation  3  Simulation  2  Simulation  1 


Figure  5.37:  Comparison  of  Reynolds  stress  data  between  spectral  element  simulations 
and  PIV  data  on  vertical  plane  Z  =  0.5. 
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Figure  5.38:  First  four  POD  modes. 
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Results  and  analysis 

In  this  section,  we  present  the  results  of  the  simulation.  We  define  the  coordinates  xs  =  x  —  11.0,  ys  = 
y  —  1.25  and  =  z.  We  use  data  along  the  slices  Syi  (y  =  1.25,  ys  =  0.0,  \  =  1.25),  SZl  (z  =  0.5,  = 

0.5,  |  =  0.25)  and  Sz 2  (z  =  1.7,  z8  =  1.7,  =  0.85)  for  our  analysis.  All  statistics  presented  in  this  section 

are  computed  based  on  80  seconds  of  data  sampled  at  50  Hz. 

Figure  5.40  plots  the  mean  velocity  profiles  along  the  midspan  slice  Syi  for  all  four  streets.  The  contour 
for  U  shows  that  the  flow  separates  from  the  top  of  the  first  block  and  seems  to  reattach  to  the  top  of  the 
second  block.  The  region  >1.3  seems  to  be  relatively  unperturbed  by  the  presence  of  the  array.  A 
strong  stream  wise  velocity  gradient  exits  in  the  shear  layer  region  0.8  <  <  1.3.  The  magnitude  of  U 

is  significantly  lower  within  the  street  and  includes  regions  of  reverse  flow.  Since  the  flow  is  symmetric 
along  the  mid-span  plane,  the  spanwise  component  of  velocity  is  negligible  in  this  plane.  The  contour  for 
W  shows  a  strong  downward  draft  on  the  windward  side  of  each  street  and  an  upward  draft  on  the  leeward 
side  of  each  street.  These  contours  indicate  the  presence  of  a  single  large  recirculation  region  for  streets  2 
through  4  which  is  expected  for  the  skimming  flow  regime.  A  notable  exception  is  street  1  which  shows  two 
large  counter-rotating  recirculation  zones. 

Figures  5.41  and  5.42  plot  the  mean  velocity  profiles  along  slices  SZl  and  Sz 2  respectively.  In  both 
cases,  the  magnitude  of  stream  wise  velocity  is  greatest  at  ^  =  0  and  ^  =  2.5  respectively.  We  can  also 
notice  that  the  magnitude  of  this  velocity  decreases  consistently  as  we  move  from  street  1  to  street  4.  A 
region  of  reverse  flow  exists  near  the  center  of  each  street  on  slice  SZl .  The  contour  of  V  on  slice  SZl  shows 
fluid  being  ejected  from  the  street  near  the  edges  of  the  upstream  block.  This  pattern  is  consistent  with  the 
structure  of  the  leg  of  a  an  arch  vortex.  The  contours  for  W  on  slices  SZ2  indicate  a  strong  downdraft  in 
the  windward  region  of  each  street  and  an  updraft  in  the  leeward  region  of  each  street  consistent  with  earlier 
observations.  Also  noticeable  is  the  decrease  in  the  magnitude  of  the  downward  and  upward  drafts  as  we 
move  from  street  1  to  street  4.  Similar  patterns  are  observed  for  SZl  for  streets  2  through  4.  However  the  W 
velocity  contour  on  slice  SZ2  seems  to  indicate  a  second  counter-rotating  recirculation  zone  close  to  the  wall 
consistent  with  the  earlier  observations. 
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Figure  5.40:  Mean  velocity  and  pressure  contours  on  slice  Syi. 
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Figure  5.41:  Mean  velocity  and  pressure  contours  on  slice  SZl . 


(d) 


Figure  5.42:  Mean  velocity  and  pressure  contours  on  slice  SZ2. 


Figure  5.43  plots  the  iso-contours  of  Ti  (Eq.  4.2)  indicating  the  locations  of  cores  of  arch  vortices  in  the 
urban  street  canyon  for  0°  angle  of  incidence.  A  threshold  value  of  0.4  is  used  for  generating  these  plots.  We 
can  see  that  for  zero  angle  of  incidence,  arch  vortices  are  located  symmetrically  with  respect  to  slice  Syi . 
We  can  also  notice  the  horse- shoe  vortex  in  front  of  the  first  block  and  vortices  on  top  and  sides  of  the  first 
block. 
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Figure  5.43:  Isocontours  of  Ti  as  defined  in  Chapter  4  (Eq.  4.2)  indicating  the  locations  of  cores  of  arch 
vortices  in  the  urban  street  canyon  for  0°  angle  of  incidence. 


Figure  5.44  plots  the  urms ,  vrms ,  wrms  and  turbulent  kinetic  energy  contours  on  slice  Syi  for  all  four 
streets.  We  can  see  that  the  intensity  of  turbulent  kinetic  energy  is  high  in  the  roof  shear  layer  region  and 
in  the  windward  side  of  each  street  just  upstream  of  the  downstream  block.  There  is  a  steady  decrease 
in  the  magnitude  of  the  turbulent  kinetic  energy  as  we  move  from  street  1  to  street  4.  In  the  shear  layer, 
the  primary  contribution  to  the  turbulent  kinetic  energy  comes  from  the  streamwise  component(^rms)  of 
velocity  fluctuations  with  a  smaller  contribution  from  the  wall  normal  component  (wrms).  On  the  windward 
side  of  the  street,  the  span  wise  component  of  velocity  fluctuations  (vrms)  is  the  major  contributor  to  the 
turbulent  kinetic  energy  with  a  minor  contribution  from  the  wall  normal  component  ( wrrns )  of  fluctuating 
velocity. 

Figures  5.45  and  5.46  plot  the  contours  of  urms ,  vrmS9  wrms  and  turbulent  kinetic  energy  on  slices  SZl 
and  Sz 2  respectively.  In  this  case,  the  intensity  of  turbulent  kinetic  energy  is  high  in  the  side  shear  layers  and 
on  the  windward  side  of  each  street.  The  intensity  of  turbulent  kinetic  energy  decreases  rapidly  as  we  move 
from  street  1  to  street  4.  The  streamwise  component  of  velocity  fluctuations  ( urrns )  is  the  major  contributor 
to  turbulent  kinetic  energy  in  the  side  shear  layers  whereas  the  spanwise  component  of  velocity  fluctuations 
( Vrms )  is  the  biggest  contributor  to  turbulent  kinetic  energy  in  the  windward  side  of  of  each  street. 

Finally,  Figures  5.47,  5.48  and  5.49  plot  contours  of  Reynolds  stresses  on  slices  Syi,  SZl  and  Sz 2  re¬ 
spectively.  On  slice  Syi  the  <  uw  >  component  of  the  Reynolds  stress  tensor  is  dominant  in  the  top  shear 
layer.  The  Reynolds  stresses  are  small  in  other  regions.  On  slices  SZl  and  SZ2,  the  <  uv  >  component  of 
the  Reynolds  stress  tensor  is  dominant  in  the  side  shear  layers.  The  strength  of  the  shear  layer  decreases 
steadily  as  we  move  from  street  1  to  street  4. 
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Figure  5.44:  urms ,  vrms ,  wrms  and  turbulent  kinetic  energy  contours  on  slice  Syi . 


106 


Figure  5.45:  urrns ,  vrrns ,  wrrns  and  turbulent  kinetic  energy  contours  on  slice  SZl . 


Figure  5.46:  urrns ,  vrrns ,  wrrns  and  turbulent  kinetic  energy  contours  on  slice  SZ2. 
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Figure  5.47:  Contours  of  (uv),  (uw)  and  (vw)  on  slice  Syi. 
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5.2.3  Numerical  simulation  of  the  urban  boundary  layer 
at  15°  angle  of  incidence 

For  this  simulation,  the  domain  consists  of  the  full  7x5  array  of  blocks  with  H/W  =  2,  S/W  =  1.5.  The 
flow  is  incident  at  an  angle  of  15°  to  the  array.  The  Reynolds  number  ( Ren ),  based  on  block  height  and  inlet 
velocity  at  block  height,  is  Ren  =  UhH/v  =  6283.  The  inflow  boundary  is  located  approximately  10  block 
widths  upstream  of  the  first  block  and  the  outflow  boundary  is  located  40  block  widths  downstream  of  the 
last  block.  A  symmetry  boundary  is  condition  is  used  on  the  top  and  is  located  at  17  block  widths  away  from 
the  bottom  wall.  A  wall  boundary  condition  is  used  in  the  span  wise  direction.  These  walls  are  about  16  block 
widths  away  from  the  outermost  block.  Inlet  velocity  profiles  obtained  from  the  hotwire  data  provide  the 
inlet  boundary  conditions  for  the  simulation.  152936  spectral  elements  with  an  eighth-order  approximation 
within  each  spectral  element  are  used  for  this  simulation.  This  corresponds  to  about  111  million  degrees  of 
freedom.  The  backward  difference  scheme  is  used  for  diffusion  terms  and  extrapolation  is  used  for  nonlinear 
terms.  De-aliasing  and  low  pass  filtering  is  used  for  stabilizing  the  simulation.  This  simulation  is  run  on 
the  MIRA  supercomputer  at  ALCF.  8192  cores  with  two  threads  per  core  are  used  for  this  simulation.  This 
corresponds  to  16384  processors.  The  code  scales  to  about  10,000  degrees  of  freedom  per  processor.  A 
cross-section  of  the  mesh  used  for  this  simulation  is  shown  in  figure  5.50.  This  simulation  is  still  in  progress 
and  the  statistics  presented  here  are  based  on  6.2  seconds  sampled  at  20  Hz.  While  this  sample  size  is  not 
sufficient  to  compute  good  quality  statistics,  we  can  see  some  of  the  basic  changes  in  flow  structure  for 
non-zero  angle  of  incidence. 

Figure  5.51  plots  the  mean  velocity  profiles  along  the  midspan  slice  Syi  for  all  four  streets.  The  contour 
for  U  shows  that  the  size  of  the  shear  layer  region  is  larger  for  this  case  compared  to  0°  angle  of  incidence. 
The  contour  for  span  wise  velocity  V  is  significantly  different  from  the  0° -angle  of  incidence  case.  We  can 
see  a  strong  channeling  effect  in  the  windward  region  of  each  street.  The  flow  is  from  the  positive  y- axis  to 
the  negative  y- axis.  The  contour  for  W  shows  a  downward  draft  on  the  windward  side  of  each  street  and  an 
upward  draft  on  the  leeward  side  of  each  street.  This  is  similar  to  the  0° -angle  of  incidence  case. 
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Figure  5.50:  Cross-section  of  mesh  along  slice  SZl  . 
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Figure  5.51:  Mean  velocity  and  pressure  contours  on  slice  Syi  for  15°  angle  of  incidence. 


Figures  5.52  and  5.53  plot  the  mean  velocity  profiles  along  slices  SZl  and  Sz 2  respectively.  In  both 
cases,  the  stream  wise  velocity  distribution  loses  symmetry  with  respect  to  slice  Syi .  The  span  wise  velocity 
component  V  is  strong  in  the  windward  region  of  the  street  and  its  magnitude  reduces  from  street  1  to 
street  4.  This  loss  of  symmetry  is  also  noticed  for  the  W  component. 
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Figure  5.52:  Mean  velocity  and  pressure  contours  on  slice  SZl  for  15°  angle  of  incidence. 


Figure  5.53:  Mean  velocity  and  pressure  contours  on  slice  Sz 2  for  15°  angle  of  incidence. 
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Figure  5.54:  Isocontours  of  Ti  as  defined  in  Chapter  4  (Eq.  4.2)  indicating  the  locations  of  cores  of  arch 
vortices  in  the  urban  street  canyon  for  15°  angle  of  incidence. 


Chapter  6 

Discussion 


6.1  Arch  vortex 

6.1.1  Comparison  experimental  and  numerical  results 

Figure  6.1  presents  a  side-by-side  comparison  of  the  arch  vortices  represented  by  iso-surfaces  of  Ti  from 
the  experimental  results  presented  in  Chapter  4  and  from  the  numerical  results  presented  in  Chapter  5  for  the 
AOI  of  0° .  In  both  cases  the  arch  vortices  are  well  captured  and  the  best  agreement  between  experimental 
and  numerical  results  is  observed  in  streets  3  and  4. 

6.1.2  Arch  vortex  core  location  dependence  on  street  and  angle  of  incidence 

Three-dimensional  T  i  iso- surfaces  presented  earlier  provide  general  information  about  the  location  of  the 
vortical  structures;  however  to  more  precisely  present  comparisons  regarding  the  location  of  the  arch  vortex 
legs,  we  focus  on  the  maximum  of  Ti  which  represents  the  actual  center  of  vortical  structures.  To  illustrate 
this,  T i  contours  are  presented  in  Figure  6.2  with  the  velocity  vector  field  superimposed  for  the  AOI  of  0° 
case  in  the  first  street.  The  data  are  presented  in  an  Sz  plane  at  z/H  =  0.5.  The  maxima  of  Ti  are  shown  as 
magenta  stars  on  top  of  the  T i  contours.  The  detection  of  the  maximum  T i  corresponds  well  with  the  center 
of  the  vortices  shown  with  the  2D  vector  field  in  the  displayed  Sz  plane.  The  accuracy  of  this  core  detection 
method  which  is  defined  based  on  the  data  acquisition  grid  dimensions,  is  1  mm  in  the  x  direction  and  2.5 
mm  in  span  ( y ). 

In  Figure  6.3,  we  use  the  same  vortex  core  location  detection  technique  but  we  now  present  a  contour  of 
the  TKE  instead  of  Ti  in  order  to  illustrate  the  spatial  connection  between  the  arch  vortex  and  the  regions 
dominated  by  the  turbulence.  Figure  6.3  presents  the  vortex  core  locations  for  all  four  streets  and  the  four 
AOI  cases  investigated.  We  also  extract  the  actual  angle  of  the  arch  vortex  structure  with  respect  to  an  axis 
aligned  with  the  7/-axis,  depicted  by  the  magenta  lines.  A  summary  of  these  angles  is  given  in  Table  6.1. 
Figure  6.3(a)  shows  that  for  AOI  =  0°,  the  arch  vortex  moves  downstream  from  street  1  to  street  4.  For  AOI 
=  0°  the  arch  vortex  angle  with  respect  to  the  street  should  be  equal  to  zero  in  each  street.  The  variation  in 
the  measured  angles  for  all  four  streets,  listed  in  Table  6.1,  arises  from  the  slight  asymmetry  in  the  incoming 
flow  triggering  a  channeling  effect  which  is  known  to  be  strong  even  for  small  AOIs  (see  Monnier  et  al. 
[2010]).  For  15°  AOI,  see  Figure  6.3(b),  the  arch  vortices  are  significantly  tilted  within  the  streets.  Both 
legs  are  still  within  the  region  resolved  with  the  SPIV  measurements.  It  can  be  seen  that  one  leg  is  leaving 
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(a)  Experimental 


(b)  Numerical 

Figure  6.1:  Arch  vortex,  using  Ti  Iso-surfaces  for  AOI  =  0°  for  both  experimental  and  numerical  results. 
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Figure  6.2:  Mean  velocity  vector  field  (black  vectors),  Ti  magnitude  contour  background,  vortex  core  loca¬ 
tions  (magenta  stars),  xy  slice  at  z/H  =  0.5. 


the  region  directly  behind  the  first  block  in  street  1  (—0.5  <  y/w  <  0.5)  and  entering  the  intersection.  For 
streets  2  to  4,  both  legs  are  still  lying  in  the  —0.5  <  y/w  <  0.5  region.  As  the  angle  is  increased  to  30°, 
Figure  6.3(c),  a  single  leg  is  captured  within  streets  1  and  2.  The  second  leg,  if  it  exists,  would  most  likely 
have  moved  out  of  the  region  resolved  in  our  measurements.  This  is  a  point  that  we  will  be  able  to  address 
in  the  near  future  with  a  numerical  simulation  of  this  case.  Within  streets  3  and  4,  both  legs  are  captured 
within  the— 0.5  <  y/w  <  0.5  region.  The  angle  of  the  arch  vortex  with  respect  to  the  street  is  larger  than 
for  the  15°  AOI  configuration,  as  can  be  seen  in  Table  6.1.  Finally,  for  the  45°  AOI,  a  single  leg  is  captured 
in  each  street.  Again,  our  assumption  is  that  a  second  leg  would  exist  in  the  region  not  resolved  by  our 
measurements. 

A  comparison  between  experimental  and  numerical  results  is  presented  in  Figure  6.4.  The  agreement 
between  the  two  is  fairly  good  in  streets  2,  3  and  4.  Differences  in  the  location  of  the  arch  vortex  are 
observed  in  street  1  but  the  TKE  distribution  is  very  similar  in  both  data  sets.  Apart  from  street  1,  the  angle 
of  the  arch  vortex  with  respect  to  the  street  obtained  from  the  numerical  simulation  is  very  comparable  to 
the  experimental  results,  see  Table  6.1. 

Table  6.1:  Arch  horizontal  axis  inclination  (tp)  in  degree  at  z/H=0.5. 


Street  # 

1st 

2nd 

3rd 

4th 

Numerical 

AOI=0° 

20.0° 

-3.5° 

-2.9° 

9.5° 

Experimental 

AOI=0° 

4.5° 

-10.2° 

2.6° 

b 

o 

AOI=15° 

29.6° 

39.7° 

37.9° 

29.8° 

AOI=30° 

- 

- 

50.2° 

51.9° 

AOI=45° 

- 

- 

- 

- 
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(a)  AOI  =  0° 


(b)  AOI  =  15° 


(c)  AOI  =  30°  (d)  AOI  =  45° 


Figure  6.3:  Mean  velocity  vector  field  (black  vectors),  TKE  magnitude  contour  background,  vortex  core 
locations  (magenta  stars),  xy  slice  at  z/H  =  0.5. 


(b)  AOI  =  0°,  numerical 


Figure  6.4:  Mean  velocity  vector  field  (black  vectors),  TKE  magnitude  contour  background,  vortex  core 
locations  (magenta  stars),  xy  slice  at  z/H  =  0.5  for  both  experimental  and  numerical  results  for  AOI  =  0°. 
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6.2  Turbulence  characteristics  of  the  gusts  within  the  streets 

In  this  section,  we  are  focusing  our  attention  on  the  probability  density  functions  (PDFs)  associated  with  the 
fluctuating  part  of  the  velocity  field  (u',vf  and  w').  From  the  experimental  data,  we  extract  the  instantaneous 
velocity  components,  u ',  v'  and  w'  from  the  SPIV  snapshots  at  two  specific  positions  in  space  which  are 
chosen  for  their  large  contribution  to  the  overall  TKE.  The  first  position  is  depicted  in  Figure  6.5(a)  as 
a  black  square  in  an  xy  slice  at  z/H  =  0.25  and  its  coordinates  are  xj S  =  0.5,  y/W  =  —0.5  and 
z/H  =  0.25.  Similarly,  the  second  position  is  shown  in  Figure  6.5(b)  in  an  x-z  slice  with  coordinates 
equal  to  x/S  =  0.5,  y/W  =  0,  z/H  =  1.  These  two  locations  are  closely  connected  with  the  shear  layers 
forming  off  the  top  and  sides  of  the  buildings  where  significant  contributions  to  the  TKE  are  observed.  Also 
included  in  the  figures  are  probability  density  functions  of  the  gusts;  that  is,  the  PDFs  of  the  individual 
fluctuating  velocity  components  (normalized  by  Uh)  at  the  position  in  space  indicated  by  the  black  square. 
Since  the  choice  of  the  black  square  was  based  on  the  relatively  high  level  of  TKE,  we  expect  large  gust 
amplitudes  in  the  corresponding  PDFs.  It  is  worth  noting  here  that  instantaneous  gust  magnitudes  can  reach 
reach  between  20%  to  40%  of  the  incoming  wind  velocity  as  normalized  by  the  roof  level  Uh  as  shown  in 
the  plots.  Also  included  in  these  plots  are  the  three  correlation  coefficients,  Ruv ,  Ruw  and  Rvs  which  give 
an  indication  of  the  coupling  between  the  different  gust  components. 

Figure  6.6  presents  a  comparison  of  the  gust  PDFs  in  street  3  for  the  four  AOIs  investigated.  In  addition  to 
the  PDFs,  a  Gaussian  curve  was  fitted  to  each  PDF  and  the  corresponding  mean,  /i,  and  standard  deviation,  a 
for  each  fluctuating  velocity  component  are  shown  along  with  the  three  correlation  coefficients.  Comparing 
the  PDFs  of  uf,  v'  and  w'  for  all  four  AOIs,  it  can  be  seen  that  the  overall  distribution  of  the  gust  is  not 
noticeably  affected  by  the  wind  direction.  The  only  significant  difference  is  observed  in  street  4  where 
the  correlation  coefficient  Ruv  drops  from  about  0.6  to  0.2  indicating  a  decoupling  of  the  streamwise  and 
span  wise  gusts.  The  overall  similarity  of  the  gust  profiles  in  the  x-y  slice  at  z/H  =  0.25  is  most  likely  due 
to  the  fact  that  this  area  of  the  street  is  more  shielded  from  the  incoming  wind. 

When  looking  at  the  second  location,  at  roof  level  of  the  blocks,  the  effect  of  the  AOI  is  more  evident  as 
shown  by  the  results  given  in  Figure  6.7.  At  an  AOI  of  0°,  the  u'  component  PDF  has  a  larger  cru>  =  0.20  as 
compared  with  the  v'  or  w'  distributions.  As  the  AOI  is  increased,  the  spanwise  v'  component  distribution 
gets  wider  and  eventually  becomes  larger  than  its  u'  counterpart.  In  addition,  the  PDFs  are  getting  more 
skewed  as  AOI  is  increased.  The  roof-level  region  is  much  more  sensitive  to  the  wind  direction  as  compared 
with  the  location  closer  to  the  ground  discussed  above.  In  terms  of  correlation  coefficients  between  any 
two  gust  components,  the  Ruw  is  dominant  for  the  0°  AOI  case.  As  the  AOI  is  increased,  the  streamwise 
gust  gains  correlation  with  the  other  two  gust  components.  By  45°,  the  dominant  correlation  coefficient  is 
Ruv  but  both  Ruw  and  Rvw  are  also  significantly  larger.  The  wind  direction  has  the  effect  of  redirecting 
the  gust  and  redistributing  the  gusts  in  all  three  directions  with  a  significant  coupling  between  the  different 
components. 

We  now  perform  a  side-by-side  comparison  of  PDFs  between  the  experimental  data  and  the  numerical 
data.  Figure  6.8  presents  such  a  comparison  for  the  first  two  streets  at  the  x/S  =  0.5,  y/W  =  —0.5, 
z/H  =  0.25  location  (near  the  ground)  for  the  0°  AOI  while  Figure  6.9  presents  the  same  comparison  for 
streets  3  and  4.  Apart  from  some  differences  in  street  1,  both  the  standard  deviations  of  the  fluctuating 
velocity  components  and  the  correlation  coefficients  are  very  close  between  the  experimental  and  numerical 
results,  indicating  that  the  turbulent  characteristics  of  the  flow  field  are  well  captured  by  the  numerical 
simulation.  This  last  point  is  also  true  when  comparing  the  gust  PDFs  at  the  second  location  (at  roof  level) 
as  is  illustrated  by  Figures  6.10  and  6.11. 
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(a)  xy  slice  at  z/H  =  0.25  showing  urms  contours  and  a  sample  of  gust  probability  density  function  based  on  data 
extracted  from  a  point  at  (x/S  =  0.5,  y/W  =  —0.5,  z/H  =  0.25) 


(b)  xz  slice  at  y/W  =  0  showing  urms  contours  and  a  sample  of  gust  probability  density  function  based  on  data 
extracted  from  a  point  at  (x/S  =  0.5,  y/W  =  0,  z/H  =  1) 

Figure  6.5:  Spatial  locations  (black  squares)  used  to  extract  gusts  probability  density  functions 
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(a)  AOI  =  0°,  experimental  results  in  street  3  (b)  AOI  =15°,  experimental  results  in  street  3 


(c)  AOI  =  30°,  experimental  results  in  street  3  (d)  AOI  =  45°,  experimental  results  in  street  3 


Figure  6.6:  PDFs  of  gusts  at  a  specific  point  in  the  street:  (x/S  =  0.5,  y/W  =  —0.5,  z/H  =  0.25)  for  all  4 
AOIs. 
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percentage  of  time  percentage  of  time 


(a)  AOI  =  0°,  experimental  results  in  street  3  (b)  AOI  =15°,  experimental  results  in  street  3 


(c)  AOI  =  30°,  experimental  results  in  street  3  (d)  AOI  =  45°,  experimental  results  in  street  3 


Figure  6.7:  PDFs  of  gusts  at  a  specific  point  in  the  street:  (x/S  =  0.5,  y/W  =  0,  z/H  =  1)  for  all  4  AOIs. 
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(a)  AOI  =  0° ,  experimental  results  in  street  1 


(b)  AOI  =  0° ,  numerical  results  in  street  1 


(c)  AOI  =  0° ,  experimental  results  in  street  2 


(d)  AOI  =  0° ,  numerical  results  in  street  2 


Figure  6.8:  PDFs  of  gusts  at  a  specific  point  in  the  street:  (x/S  =  0.5,  y/W  =  —0.5,  z/H  =  0.25) 
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percentage  of  time  percentage  of  time 


(a)  AOI  =  0° ,  experimental  results  in  street  3 


(b)  AOI  =  0° ,  numerical  results  in  street  3 


(c)  AOI  =  0° ,  experimental  results  in  street  4  (d)  AOI  =  0° ,  numerical  results  in  street  4 

Figure  6.9:  PDFs  of  gusts  at  a  specific  point  in  the  street:  (x/S  =  0.5,  y/W  =  —0.5,  z/H  =  0.25) 
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percentage  of  time  percentage  of  time 


(a)  AOI  =  0° ,  experimental  results  in  street  1  (b)  AOI  =  0°,  numerical  results  in  street  1 


(c)  AOI  =  0° ,  experimental  results  in  street  2  (d)  AOI  =  0° ,  numerical  results  in  street  2 


Figure  6.10:  PDFs  of  gusts  at  a  specific  point  in  the  street:  (x/S  =  0.5 ,y/W  =  0,  z/H  =  1) 
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percentage  of  time  percentage  of  time 


(a)  AOI  =  0° ,  experimental  results  in  street  3 


(b)  AOI  =  0° ,  numerical  results  in  street  3 


(c)  AOI  =  0° ,  experimental  results  in  street  4 


(d)  AOI  =  0° ,  numerical  results  in  street  4 


Figure  6.11:  PDFs  of  gusts  at  a  specific  point  in  the  street:  (x/S  =  0.5 ,y/W  =  0,  z/H  =  1) 
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6.3  Connection  between  arch  vortices  and  high  turbulence  areas 

In  this  section,  we  are  looking  at  the  arch  vortex  location  with  respect  to  regions  of  high  turbulence  (large 
Urms ,  Vrms  or  wrrns).  To  do  so,  we  are  presenting  3D  iso-surfaces  of  Ti,  urrns  and  vrrns  in  Figure  6.12.  The 
experimental  data  are  presented  in  street  1  for  the  four  AOIs  investigated.  The  thresholds  used  to  plot  the 
regions  of  high  turbulence  are  selected  as  75%  of  the  local  maximum  of  urrns  and  vrms  within  the  street.  For 
the  0°  AOI,  see  Figure  6.12(a),  the  two  regions  of  large  urrns  are  essentially  due  to  the  shear  layers  formed 
off  the  sides  of  the  blocks  with  the  arch  vortex  being  “trapped”  between  these  two  regions.  Similarly,  the 
region  of  large  vrrns  observed  near  the  windward  face  of  the  downstream  block  is  located  downstream  of  the 
arch  vortex.  The  large  spanwise  fluctuations  associated  with  vrms  are  observed  in  the  middle  of  the  street 
in  the  spanwise  direction  (y/W  =  0)  which  is  also  between  the  legs  of  the  arch  vortex.  The  arch  vortex  is 
surrounded  by  regions  of  high  turbulence  but  its  core  sits  in  a  low  turbulence  area.  As  the  AOI  is  increased  to 
15°,  see  Figure  6.12(b),  the  arch  vortex  is  tilted  with  respect  to  the  street  axis.  The  regions  of  high  turbulence 
are  also  redistributed  around  the  arch  vortex.  The  spanwise  turbulence  region  is  shifted  and  again  aligned 
with  the  middle  of  the  two  legs  of  the  arch  vortex.  For  the  30°  AOI,  see  Figure  6.12(c),  the  picture  is  very 
similar  to  the  15°  case.  In  the  45°  AOI  case,  the  trend  is  similar  but  an  additional  region  of  high  spanwise 
turbulence  appears  due  to  the  side  shear  layer  with  the  arch  vortex  still  sitting  in  the  low  turbulence  region. 

Comparing  the  experimental  results  with  the  numerical  results,  see  Figure  6.13,  the  same  trends  are 
observed.  The  main  difference  lies  in  the  shape  of  the  arch  vortex  which  partly  explains  the  difference 
observed  in  the  fluctuating  velocity  PDFs  presented  in  the  previous  section. 
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(a)  AOI  =  0° 


(b)  AOI  =  15° 


75%  Of 

75%  of 


(c)  AOI  =  30° 


(d)  AOI  =  45° 


Figure  6.12:  Arch  vortex  and  regions  of  high  urrns  and  vrms. 


132 


Figure  6.13:  Arch  vortex  and  regions  of  high  urrns  and  vrrns  for  both  experimental  and  numerical  results. 
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Chapter  7 


Conclusions 


7.1  Experiment 

In  this  study,  a  scaled  urban  environment  model  is  simulated  in  a  low-speed  wind  tunnel.  Stereoscopic 
particle  image  velocimetry  (SPIV)  was  used  to  provide  three-dimensional  velocity  data  in  the  urban  streets 
to  achieve  a  thorough  spatial  map  of  the  mean  and  turbulent  flow.  Contour  slices  and  iso-surfaces  are  the 
tools  used  to  document  the  complex  flow  in  four  consecutive  streets  within  the  urban  array.  Additionally, 
the  effect  of  the  incoming  wind  direction  on  the  flow  characteristics  was  studied  for  incidence  angles  of  0° 
to  45°. 

The  channelling  of  the  flow  down  the  streets  is  easily  seen  by  the  results  for  the  non-zero  incidence 
angles  of  15°,  30°  and  45°.  The  channelling  of  the  flow  is  also  reflected  in  the  tilting  of  the  arch  vortex.  For 
the  0°  case  the  arch  vortex  in  street  1  is  located  relatively  close  to  the  upstream  building  and  its  vertical  axis 
is  almost  perpendicular  to  the  ground.  Further  downstream  in  the  array  the  arch  vortex  is  observed  to  tilt 
(around  the  spanwise  coordinate)  in  the  downstream  direction.  T i  was  used  to  locate  the  core  locations  of 
the  arch  vortex  legs  in  the  horizontal  streamwise-spanwise  plane  at  a  wall-normal  position  of  z/H  =  0.5. 
These  results  show  that  there  is  also  an  effect  of  AOI  on  the  tilting  of  the  arch  vortex  around  the  streamwise 
coordinate. 

The  TKE  results  for  all  four  AOI  cases  investigated  illustrate  that  the  flow  through  the  array  seems  to 
reach  an  equilibrium  condition  in  as  little  as  3  to  4  streets.  The  turbulence  levels  from  street  1  to  2  were 
shown  to  decrease  significantly  with  smaller  decreases  from  street  to  street  after  that.  Regarding  the  ef¬ 
fect  of  AOI  on  the  turbulence  levels  it  is  observed  that  the  TKE  level  increases  slightly  from  AOI=0°  to 
AOI=15°  and  then  decreases  as  the  AOI  increases  to  30°  and  45°.  By  decomposing  the  TKE  into  its  three 
components  it  is  shown  that  the  streamwise  ( u /2)  and  spanwise  ( v /2)  fluctuating  components  are  more  signif¬ 
icant  than  the  wall-normal  ( w /2)  component.  The  relationship  between  the  turbulence  and  gust  magnitudes 
was  investigated  by  looking  at  the  probability  density  functions  for  all  three  velocity  components  in  regions 
corresponding  to  high  TKE.  Gust  magnitudes  as  large  as  20-40%  of  Uh  in  all  three  velocity  components 
are  common  for  all  cases  considered.  The  correlations  between  the  gust  components  were  also  studied  and 
it  was  found  that  lower  in  the  street  the  correlations  between  the  streamwise  and  spanwise  gusts,  (RUv)> 
were  dominant  for  all  four  AOI  cases  and  the  two  other  cross  correlations  (spanwise/wall-normal  ( Rvw  )  & 
streamwise/wall-normal  (Ruw))  were  negligible.  At  higher  wall-normal  positions  in  the  array  (z/H  =  1) 
the  Rvw  &  Ruw  correlations  increased  as  compared  with  the  lower  z/H  =  0.25  condition.  At  z/H  =  1  for 
AOI  =  0°  the  Ruv  correlation  is  negligible  whereas  the  correlation  between  the  streamwise  and  wall-normal 
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gusts  ( Ruw )  is  relatively  large.  As  AOI  increases  at  this  z/H  =  1  position  Ruw  decreases  as  AOI  increases 
and  for  AOI  values  of  30°  and  45°  all  three  correlations  are  significant. 

Overlaying  the  arch  vortex  structures  with  the  TKE  components  reveals  that  the  vortex  cores  are  located 
in  regions  of  relatively  low  turbulence  but  surrounded  by  high-turbulence  regions.  For  AOI  =  0°  the  high 
urrns  regions  are  due  to  the  shear  layers  from  the  sides  of  the  building.  Similarly,  the  region  of  large  vrrns 
observed  near  the  windward  face  of  the  downstream  block  is  located  downstream  of  the  arch  vortex.  As  the 
incidence  angle  is  increased  the  arch  vortex  is  tilted  with  respect  to  the  street  axis  with  the  regions  of  high 
turbulence  being  redistributed  around  the  arch  vortex.  In  addition  the  region  of  large  v rrns  values  is  shifted 
and  again  aligned  with  the  middle  of  the  two  legs  of  the  arch  vortex. 


7.2  Numerical 

As  part  of  this  project,  a  high-order  accurate  incompressible  Navier-Stokes  solver  capable  of  performing 
high  fidelity  simulations  of  flows  of  engineering  interest  was  built.  The  physical  problem  to  be  tackled  here 
required  a  solver  which  is  highly  scalable,  supports  complex  geometries  and  has  low  dispersion  and  diffusion 
errors. 

We  therefore  developed  Spec  solve,  a  scalable  spectral  element  solver.  The  solver  is  coded  in  C++  and 
uses  the  MPI  (Message  Passing  Interface)  library  for  communication.  It  exploits  the  object-oriented  features 
of  C++  and  uses  dynamic  memory  allocation  for  optimal  memory  usage.  The  code  is  mostly  self-contained 
and  the  only  external  dependencies  are  the  standard  LAPACK  and  BLAS  libraries.  It  supports  meshes 
generated  by  general-purpose  mesh  generators  like  CUBIT  and  GAMBIT.  An  efficient  mesh  partitioner, 
based  on  parallel  recursive  spectral  bisection  algorithm,  is  built  for  generating  high  quality  mesh  partitions. 

The  solver  uses  the  Vn  —  Vn-2  formulation  for  spatial  discretization  of  velocity  and  pressure.  Velocity 
is  represented  within  each  spectral  element  using  a  tensor  product  of  Lagrangian  basis  functions  based  on 
GLL  (Gauss-Lobatto-Legendre)  nodes  and  C°  continuity  is  enforced  between  adjacent  spectral  elements. 
Pressure  within  each  spectral  element  is  represented  using  a  tensor  product  of  one-dimensional  Legendre 
polynomials  of  appropriate  order  and  inter-element  continuity  is  not  explicitly  enforced.  A  semi-implicit 
scheme,  which  treats  the  Stokes  operator  implicitly  and  the  nonlinear  term  explicitly,  is  used  for  temporal 
integration.  A  fractional  step  scheme  with  second-order  temporal  accuracy  is  used  to  decouple  velocity  and 
pressure.  This  requires  the  solution  of  a  Helmholtz  system  for  each  component  of  velocity  and  consistent 
Poisson  equation  for  pressure  at  each  time  step.  While  the  Helmholtz  system  can  be  solved  efficiently  using 
the  Jacobi  preconditioned  conjugate  gradient  method,  the  solution  of  the  consistent  Poisson  equation  remains 
the  principal  bottleneck  in  the  numerical  solution  of  the  Navier-Stokes  equations. 

To  address  this,  a  multilevel  strategy  was  implemented  to  build  a  scalable  solver  for  pressure.  The 
deflation  approach  is  used  to  decompose  pressure  into  fine  and  coarse  components.  An  FDM-based  block- 
Jacobi  preconditioner  is  used  for  efficient  solution  of  the  fine  pressure  problem.  A  parallel  direct  solver 
and  a  state-of-the-art  algebraic  multigrid  solver  are  implemented  to  solve  the  coarse  pressure  problem.  The 
parallel  direct  solver  is  used  for  meshes  with  less  than  105  spectral  elements  whereas  the  algebraic  AMG 
solver  is  used  for  larger  mesh  sizes.  A  multi-threaded  C++  code  was  built  for  setting  up  the  data  needed  for 
the  AMG  solver. 

The  solver  was  tested  for  accuracy  on  various  benchmark  problems.  Two-dimensional  test  cases  com¬ 
prised  Wannier  flow,  Kovasznay  flow,  vortex  shedding  from  a  circular  cylinder  at  various  Reynolds  num¬ 
bers,  and  the  highly  sensitive  Orr-Sommerfeld  problem.  Three-dimensional  test  cases  include  the  three- 
dimensional  Kovasznay  flow  and  the  flow  over  a  backward-facing  step  at  Re  =  172  and  Re  =  343.  In  all 
cases,  the  results  are  in  excellent  agreement  with  the  corresponding  exact  solutions  predicted  by  theory  and 
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experimental  data. 

Additional  issues  arise  when  high-order  solvers  are  used  outside  their  traditional  realm  of  simulating 
canonical  flows.  Typically,  most  flows  of  engineering  interest  occur  in  relatively  complex  domains  and  it  is 
fairly  common  to  have  outflow  boundaries.  In  addition,  it  is  generally  impractical  to  have  highly  resolved 
meshes  throughout  the  flow  domain.  Simulating  high  Reynolds  number  flow  generally  requires  a  stabiliza¬ 
tion  method  in  the  absence  of  LES  models,  or  exorbitantly  high  mesh  resolutions  and  an  outflow  boundary 
condition  that  can  handle  energy  influx  to  the  domain  caused  by  strong  vortices  exiting  at  the  outflow  bound¬ 
ary.  At  high  Reynolds  numbers,  our  solver  uses  filter-based  stabilization  and  supports  a  turbulent  outflow 
boundary  condition  to  stabilize  the  flow  at  the  outflow  boundary. 

Our  solver  was  used  to  simulate  flow  in  a  model  urban  street  canyon.  The  domain  consists  of  an  array 
of  blocks,  typical  of  a  modern  urban  environment,  placed  in  the  wind  tunnel.  The  Reynolds  number  based 
on  block  height  ( H )  and  inlet  velocity  at  block  height  ( Uh )  is  Ren  =  UhvH  =  6283.  The  simulation  results 
are  in  good  agreement  with  the  PIV  data  generated  from  the  experiments  and  demonstrate  the  fitness  of  the 
solver  for  production  use. 


7.3  Future  work 

Our  immediate  objectives  are  as  follows: 

•  Further  improve  the  comparison  between  experiments  and  simulation  for  the  urban  street  canyon 
simulation.  This  would  involve  using  more  accurate  inflow  boundary  conditions. 

•  Study  the  effect  of  angle  of  attack  on  flow  in  the  urban  street  canyon.  This  would  involve  simulating 
the  entire  three  dimensional  array  as  opposed  to  simulating  a  single  row  of  blocks. 

•  Extend  this  work  from  wind  tunnel  scale  to  field  scale.  This  would  require  implementation  of  robust 
LES  models  Bou-Zeid  et  al.  [2005]  which  work  with  marginally  resolved  grids  which  only  resolve  the 
upper  edge  of  the  inertial  layer. 

•  Further  improve  the  scalability  of  the  solver. 
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